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Abstract 


This  dissertation  focuses  on  increasing  the  ability  to  detect  space  objects  and  increase 
Space  Domain  Awareness  (SDA)  with  space  surveillance  sensors  through  image  processing 
and  optical  theory.  SDA  observations  are  collected  through  ground-based  radar  and  optical 
systems  as  well  as  space  based  assets.  This  research  focuses  on  a  ground-based  optical 
telescope  system,  the  Space  Surveillance  Telescope  (SST).  By  increasing  the  number  of 
detectable  Resident  Space  Objects  (RSOs)  through  image  processing,  SDA  capabilities  can 
be  expanded.  This  is  accomplished  through  addressing  two  main  degrading  factors  present 
in  typical  SDA  sensors;  spatial  undersampling  in  the  collected  data  and  noise  models  and 
assumptions  used  in  current  algorithms.  The  assigned  cost  and  a  priori  probabilities  of 
a  Bayes  Multiple  Hypothesis  Test  (MHT)  are  investigated  in  this  dissertation  to  address 
the  spatial  undersampling.  New  algorithms  are  developed  and  tested,  and  demonstrated 
improved  detection  capabilities  at  operationally  realistic  false  alarm  rates.  Additionally,  a 
new  noise  model  is  developed  which  more  accurately  represents  the  received  noise  present 
in  data  collected  with  surveillance  telescopes  under  certain  atmospheric  conditions.  These 
algorithm  have  demonstrated  probability  of  detection  improvement  of  up  to  80  percent  in 
collected  SST  data  over  the  currently  employed  detection  techniques. 
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OPTICAL  THEORY  IMPROVEMENTS  TO  SPACE  DOMAIN  AWARENESS 


I.  Introduction 


1.1  Motivation 

The  motivation  for  this  research  is  to  increase  Space  Domain  Awareness  (SDA) 
capabilities  by  proposing  a  new  detection  algorithm  that  improves  the  ability  to  detect  space 
objects.  SDA  is  a  large  field  that  encompasses  several  components  including  the  detection, 
tracking,  and  characterization  of  space  objects.  These  objects  include  satellites,  orbital 
debris,  and  Near-Earth  Asteroids  (NEA).  SDA  is  currently  achieved  through  multiple 
methods  including  ground  and  space-based  optical  systems  and  ground-based  radars.  This 
research  effort  focuses  on  ground-based  observation  telescopes,  specifically  the  Space 
Surveillance  Telescope  (SST).  The  SST  is  a  telescope  built  by  MIT  Lincoln  Labs  under 
the  direction  of  the  Defense  Advanced  Research  Projects  Agency  (DARPA).  The  mission 
of  the  SST  program  is  “to  enable  ground-based,  broad-area  search,  detection  and  tracking 
of  small  objects  in  deep  space  for  purposes  such  as  space  mission  assurance  and  asteroid 
detection  [1].”  The  SST  is  a  large  field-of-view  telescope  with  the  ability  to  quickly  scan  the 
night  sky  and  detect  and  track  objects  in  the  Earth’s  orbit,  as  well  as  deep  space  asteroids. 

There  are  multiple  stakeholders  in  SDA  data  collection.  They  include  the  Department 
of  Defense  (DoD)  and  other  US  and  foreign  government  agencies.  These  agencies  have 
published  space  policies  that  include  improving  SDA  data  collection  as  a  priority  for 
their  organization  [2-4].  It  is  important  that  space  objects  in  orbit  are  cataloged  in  order 
to  mitigate  their  potential  to  cause  severe  damage  to  space  assets  through  collisions  or 
malicious  actions.  Avoiding  space  objects  is  critical  to  maintaining  functional  space  assets 
and  retaining  a  tactical  edge  in  the  space  domain.  Information  collected  by  the  SST  and 


1 


other  SDA  telescopes  is  used  for  achieving  the  goals  outlined  in  US  space  policies  and  for 
maintaining  a  robust  and  accurate  SDA  picture.  Any  improvements  in  the  ability  of  the 
SST  to  detect  objects  in  orbit  will  help  in  satisfying  this  goal  of  improved  detection  and 
characterization  of  space  objects.  Figure  1.1  is  a  visual  representation  of  objects  located 
in  Earth’s  orbit.  It  provides  a  graphical  example  of  the  orbital  areas  that  contain  the  most 
congestion  and  concern.  95  percent  of  the  objects  displayed  in  Figure  1.1  are  orbital  debris 
[5]. 


Figure  1.1:  NASA  visual  representation  of  objects  in  Earth’s  orbit  [5]. 


In  addition  to  the  detection  of  objects  in  orbit,  detecting  NEAs  is  another  part  of  the 
SST’s  mission.  National  Aeronautics  and  Space  Administration  (NASA)  has  been  tasked 
with  detecting  90  percent  of  the  NEAs  that  pose  a  severe  threat  to  humankind  by  2020  [6]. 
These  objects  are  defined  as  asteroids  larger  than  140m  in  diameter  that  have  a  perihelion 
distance  of  less  than  1.3  Astronomical  Units  from  the  sun.  In  a  2010  National  Research 
Council  (NRC)  report  on  the  progress  towards  this  goal,  it  was  determined  that  the  NEA 
survey  is  not  on  track  to  be  completed  by  2020  [7].  One  major  challenge  identified  is 
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the  lack  of  necessary  funding.  These  conclusions  were  further  detailed  in  a  2014  NASA 
Inspector  General  report  [8].  Due  to  a  lack  of  necessary  funding,  as  stated  in  these  reports, 
there  is  a  significant  benefit  to  utilizing  existing  telescope  systems  to  improve  SDA  data 
collection  efforts.  Improving  the  ability  of  the  SST  to  detect  dimmer  and  smaller  objects 
will  improve  the  progression  towards  the  mandated  goal,  with  little  or  no  additional  cost 
or  hardware.  At  the  time  of  this  report,  NASA  believes  they  have  currently  surveyed  10 
percent  of  the  90  percent  goal,  with  many  of  those  surveyed  objects  being  larger  than  1km. 

Multiple  studies  have  investigated  the  environmental  and  population  risk  factors 
caused  by  NEAs.  In  [9],  Chapman  demonstrates  that  impacts  from  objects  as  small  as  1km, 
many  of  which  are  still  not  cataloged  as  stated  previously,  have  the  potential  to  destabilize 
the  global  ecosystem  threatening  human  civilization.  Other  studies  explore  other  potential 
asteroid  impact  factors,  including  an  ocean  impact  study  by  Morrison  et.  al.  [10],  and  bias 
corrected  impact  prediction  data  by  Stuart  and  Binzel  [11]. 

1.2  System  Description 

This  dissertation  focuses  mainly  on  the  SST.  To  investigate  the  research  questions 
proposed,  both  simulated  and  captured  SST  data  is  processed  with  the  newly  developed 
detection  algorithms.  The  SST  system  is  a  Mersenne-Schmidt  design  and  is  currently 
located  at  the  White  Sands  Missile  Range  (WSMR)  in  New  Mexico  at  approximately 
8,000ft  elevation  [12].  The  weather  and  altitude  provide  generally  favorable  conditions 
for  data  collections  on  most  nights.  Table  1.1  shows  additional  key  system  parameters  of 
the  SST.  Figure  1.2  shows  an  optical  diagram  of  SST  as  well  as  a  model  of  the  complete 
telescope  system. 

The  SST  is  currently  anticipated  to  be  moved  to  Exmouth,  Australia.  This  will 
improve  the  global  coverage  compared  to  the  currently  implemented  Ground-based 
Electro-Optical  Deep  Space  Surveillance  (GEODSS)  network.  In  addition,  a  second 
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Table  1.1:  SST  System  Parameters.  [13] 


Parameter 

Focal  Length 

Primary  Mirror  /  Obscuration 
CCD  Pixels  Size 
Total  Number  Pixels 
Center  Wavelength 


Value 

3.5m 

3.5m  /  1.75m 
15/rm  x  15/rm 
6144x4096 
500nm 


telescope  may  be  installed  at  the  WSMR  to  replace  the  original.  Any  new  algorithms  can 
potentially  be  implemented  in  the  original  telescope,  or  updated  to  the  new  data  processing 
pipeline. 


L, 


3D  Uyout 


Figure  1.2:  Model  and  optical  diagram  of  the  Mersenne-Schmidt  design  for  the  SST.  [13] 
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1.3  Research  Goal  and  Objectives 

The  goal  of  this  research  is  to  improve  the  detection  capability  of  current  ground- 
based  SDA  optical  telescopes  with  new  detection  algorithms  by  addressing  two  main 
degrading  factors:  spatial  undersampling  and  the  received  noise  assumptions.  This 
dissertation  addresses  this  research  goal  with  the  following  research  questions.  These 
research  questions  led  to  the  three  main  topics  covered  in  this  dissertation. 

1.  Will  a  new  realistic  cost  function  improve  the  detection  performance  of  a  Bayes 
Criterion  MHT? 

2.  Do  the  assignments  of  a  priori  probabilities  in  a  MHT  improve  the  detection 
performance? 

3.  Can  using  a  more  realistic  noise  model  for  detection  algorithms  increase  the  ability 
to  detect  space  objects? 

Research  questions  1  and  2  address  the  spatial  undersampling,  and  research  question 
3  addresses  the  received  noise  assumptions. 

1.4  Document  Organization 

This  document  is  organized  into  6  chapters  that  cover  the  relevant  background,  as 
well  as  each  research  question  addressed  in  section  1.3.  Chapter  2  is  a  review  of  relevant 
background  and  publications  in  this  research  area.  The  background  includes  detection 
theory  and  a  review  of  current  and  newly  proposed  algorithms.  In  addition,  a  review  of 
atmospheric  theory  and  methods  for  correcting  for  the  atmospheric  effects  present  in  SDA 
data  collections  is  presented. 

Chapter  3  explores  the  development  of  a  new  detection  algorithm.  This  algorithm 
is  based  on  an  unequal-cost  MHT,  expanding  on  previously  developed  matched  filter  and 
equal  cost  MHT  algorithms. 
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Chapter  4  examines  the  segmentation  of  a  single  pixel  into  the  decision  space  for  a 
MHT.  This  segmentation  is  used  as  a  basis  for  the  prior  probabilities  used  in  a  MHT. 
Simulated  and  collected  telescope  data  are  both  analyzed. 

Chapter  5  investigates  the  underlying  statistics  of  the  received  data.  The  goal  with  this 
portion  of  the  research  is  to  expand  upon  the  Gaussian  distribution  of  noise  assumed  in  the 
Chapter  3.  The  atmosphere  and  random  photon  arrival  times  are  investigated  to  create  an 
accurate  joint  model.  Chapter  6  covers  the  future  work  that  can  be  investigated  after  this 
dissertation. 
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II.  Background  &  Literature  Review 


This  chapter  outlines  background  research  and  discoveries  that  relate  to  the  research 
topics  discussed  in  this  dissertation,  with  two  main  topic  areas.  The  first  topic  area 
is  background  on  space  object  detection,  focusing  on  both  the  history  and  current 
methodology  used.  The  second  topic  area  is  an  investigation  of  the  properties  and  effects 
of  the  atmosphere,  and  the  methods  of  dealing  with  its  impact  on  imaging  systems. 

2.1  Background  on  Space  Object  Detection 

The  SST  is  one  of  the  latest  systems  to  collect  data  on  space  objects,  including  NEA 
detection.  Before  the  use  of  Charge-Coupled  Device  (CCD)  technology,  searches  for  NEAs 
were  conducted  by  two  main  methods.  The  first  method  involved  the  human  eye  and 
using  memory  or  note-taking  to  capture  object  locations.  Through  these  observations, 
astronomers  could  determine  differences  in  object  position  over  time  in  order  to  locate 
asteroids.  When  photograph  technology  became  available,  the  comparison  of  film  images 
allowed  for  detailed  comparison  over  time.  Film  comparison  is  more  exact,  but  still 
requires  manual  analysis,  limiting  the  quantity  of  information  that  can  be  processed  and 
the  complexity  of  the  algorithms  that  can  be  used. 

Figure  2.1  is  a  plot  of  the  number  of  NEA  objects  detected  by  year.  It  demonstrates 
that  the  vast  majority  of  these  objects  have  been  found  and  cataloged  in  the  past  20  years. 
In  addition,  a  small  percentage  of  the  cataloged  objects  are  what  NASA  classifies  as  large 
NEAs. 

As  Figure  2.1  shows,  large  advancements  in  the  discovery  of  asteroids  have  occurred 
over  the  last  20  years.  These  discoveries  can  be  tied  to  the  use  of  digital  detection 
techniques.  The  Spacewatch  program  began  developing  digital  detection  techniques  in 
the  early  1980s  [14].  This  program  involves  capturing  images  with  CCDs,  and  processing 
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Figure  2.1:  Number  of  total  NEO  and  large  NEA  discoveries  as  of  15  July  2014  [8]. 


the  data  with  moving  object  detection  algorithms.  As  awareness  and  understanding  of  the 
threats  posed  by  NEAs  and  space  debris  detection  rose,  additional  programs  and  research 
efforts  were  created  to  address  these  issues.  Figure  2.2  demonstrates  the  diameters  of 
discovered  NEA  objects.  As  the  diameter  of  the  space  object  increases  it  becomes  more 
easily  detectable  with  ground-based  telescopes,  and  a  larger  percentage  of  the  objects  are 
discoverable.  Larger  objects  occur  less  frequently  than  smaller  objects,  so  there  are  less  to 
be  discovered. 

Another  program  of  note  is  the  Lincoln  Near  Earth  Asteroid  Research  (LINEAR) 
program.  Details  about  the  LINEAR  program  are  given  by  Stokes  et  al.  in  [15].  LINEAR 
uses  a  GEODSS  telescope  to  capture  images  together  with  a  Binary  Hypothesis  Test  (BHT) 
point  detector  to  detect  NEAs.  As  described  by  Viggh  et  al.,  a  point  detector  is  used  to 
generate  a  binary  map  of  ones  and  zeros  signifying  where  objects  are  present  [16].  The 
SST  currently  uses  a  similar  point  detection  algorithm  to  LINEAR.  A  binary  map  is 
generated  in  a  similar  method  to  LINEAR,  and  then  detections  are  monitored  over  three 
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Figure  2.2:  Diameter  of  currently  discovered  and  cataloged  NEAs  [8]. 


consecutive  frames  to  make  a  final  detection  decision.  If  three  consecutive  detections  are 
observed,  there  is  considered  to  be  an  object  present.  This  method  of  detection  works  well 
if  the  majority  of  the  object’s  intensity  is  in  a  single  pixel.  A  point  detector  yields  the 
largest  possible  Signal  to  Noise  Ratio  (SNR)  for  that  pixel,  and  the  greatest  chance  for 
detection.  A  point  detector  can  be  implemented  with  relative  simplicity  and  does  not  have 
large  computational  requirements.  A  major  drawback  of  this  method  is  that  the  intensity 
from  an  object  can  be  spread  over  multiple  pixels,  reducing  the  peak  intensity  of  any  single 
pixel.  Figure  2.3  is  a  block  diagram  depicting  the  entire  detection  process  used  by  LINEAR 
[16]. 

The  detection  process  contains  three  important  elements:  data  capture  and  pre¬ 
processing,  the  detection  algorithm,  and  post-processing  and  output.  When  discussing 
detection  algorithms,  the  middle  two  blocks  are  the  processes  being  referenced.  The  pre- 
and  post-processing  portions  can  remain  the  same  regardless  of  the  detection  theory  used 
to  generate  a  binary  map. 
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Figure  2.3:  LINEAR  system  detection  block  diagram  based  on  [16].  The  SST  system 
utilizes  a  similar  approach  to  detection. 


Recently  more  emphasis  has  been  put  on  the  characterization,  modeling,  and  resolved 
imaging  of  space  objects.  Gathering  this  information  can  provide  additional  detail  on  the 
object,  but  each  collection  technique  presents  its  own  technical  challenges.  This  research 
does  not  address  this  aspect  of  SDA.  The  detections  and  tracks  generated  from  the  search 
telescopes,  including  the  SST  studied  in  this  research,  will  lead  to  follow-up  investigations 
by  other  sensors  to  gather  additional  information.  The  follow-up  sensors  may  collect  the 
necessary  characterization  data.  Additionally,  the  detection  algorithms  and  techniques 
developed  in  this  dissertation  can  be  applied  to  other  research  areas.  One  example  is  star 
tracking  for  communications  or  navigation.  The  space  objects  and  stars  both  appear  as 
point  sources,  and  can  be  processed  similarly.  Increasing  the  number  of  detected  stars  in 
an  image  by  using  advanced  detection  techniques  can  potentially  increase  the  pointing  or 
navigation  accuracy  of  star  trackers. 

There  are  three  key  factors  that  impact  the  types  of  detection  algorithms  that  are 
developed,  and  their  effectiveness  in  detecting  space  objects.  It  is  important  to  understand 
the  conditions  under  which  the  data  is  collected.  The  following  assumptions  are  made 
about  the  data: 
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•  All  objects  observed  through  the  optics  are  unresolved  point  sources.  The  apparent 
size  of  space  objects  on  the  detector  are  not  necessarily  limited  to  one  pixel  due 
to  intensity  spreading  effects  including  optical  aberrations  or  atmospheric  effects. 
Since  the  objects  are  unresolved,  there  is  no  need  for  image  reconstruction  algorithms 
which  can  be  computationally  costly. 

•  The  decisions  about  an  object  being  present  or  not  are  made  only  by  the  data 
available  in  a  single  frame  or  captured  image.  The  persistence  of  Resident  Space 
Objects  (RSOs)  may  be  noted  frame  to  frame  to  reduce  false  alarms,  as  is  the  case 
with  the  SST.  There  are  algorithms  that  can  take  advantage  of  multiple  frames  to 
increase  detection  performance,  but  these  are  not  investigated  in  this  dissertation. 

•  In  Chapters  3  and  4,  the  frames  are  collected  with  long  exposure  imaging.  According 
to  Goodman,  long  exposure  images  are  captured  with  an  integration  period  of  greater 
than  10ms  [17].  Additionally,  the  integration  time  is  not  long  enough  to  cause  objects 
moving  at  sidereal  rate  to  streak. 

Noting  these  assumptions,  more  specific  detection  algorithms  can  be  developed  that 
take  advantage  of  the  known  data  collection  methods.  With  the  increase  in  computer 
memory  and  processing  speed,  increasingly  advanced  detection  algorithms  have  become 
possible.  Methods  that  do  more  than  compare  single  pixels  against  a  threshold  provide 
greater  detection  ability  at  the  cost  of  processing  complexity.  Using  the  expected  image  of 
an  object  viewed  as  a  point  source,  or  Point  Spread  Function  (PSF),  to  search  and  make 
detection  decisions  is  known  as  a  matched  filter  technique.  This  method  allows  more  than 
a  single  pixel  to  be  used  in  the  detection  decision.  Matched  filtering  effectively  averages 
the  noise  over  all  the  pixels  used.  One  current  standard  for  a  matched  filter  algorithm  used 
in  multiple  NEA  programs  is  SExtractor,  proposed  by  Bertin  in  [18].  This  is  a  software 
package  that  processes  astronomical  images  and  performs  detection  and  classification  of 
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objects.  The  portion  of  SExtractor  that  relates  to  this  paper  is  the  detection,  or  thresholding, 
step.  The  authors  propose  a  method  of  convolution  between  the  received  data  and  the 
PSF  for  faint  unresolved  objects,  the  type  of  objects  investigated  in  this  paper.  Additional 
methods  by  Bertin  applying  SExtractor  and  improving  the  algorithm  can  be  found  in 
[19,  20],  these  improvements  utilize  the  traditional  matched  filter  technique. 

Matched  filter  space  object  detectors  can  be  separated  into  two  categories:  spatial  only, 
and  spatial  and  temporal.  In  a  spatial  only  target  detection,  only  the  spatial  characteristics 
of  the  object  being  investigated  are  utilized  in  the  detection  process,  which  can  include 
the  shape  and  intensity  distribution.  This  type  of  algorithm  is  used  when  the  object  does 
not  move  significantly  during  the  integration  time  of  the  image.  Matched  filter  algorithms 
accounting  for  both  space  and  time  are  also  utilized  for  detection  of  space  objects.  In  these 
algorithms  by  Gural  et.  al  and  Pohlig,  the  spatial  and  temporal  characteristics  are  both 
used  to  make  detection  decisions  [21,  22].  These  detection  algorithms  are  not  investigated 
further  in  this  paper  because  they  do  not  match  the  data  collection  methods  used  in  the  SST, 
where  the  integration  time  does  not  allow  for  significant  orbital  motion. 

In  spatial  matched  filtering,  it  is  important  to  have  an  accurate  model  or  prediction  of 
what  the  object  is  expected  to  look  like  in  the  imaging  system.  In  [23],  O’Dell  shows  that 
the  spatial  sampling  of  the  CCD  pixels  impacts  the  resulting  image.  The  author  investigates 
the  effects  of  sampling  at  both  Rayleigh  and  Nyquist  rates,  and  demonstrates  the  impact  on 
detection  performance  of  a  matched  filter  in  undersampled  systems. 

In  [24],  Zingarelli  demonstrates  a  method  for  increasing  the  detection  performance 
over  a  traditional  matched  filter  techniquei  in  undersampled  systems.  Recognizing  that 
where  the  object  is  formed  within  a  CCD  pixel  will  change  the  distribution  of  the  image, 
the  author  creates  multiple  potential  PSFs  for  comparison.  This  leads  to  a  MHT,  where 
multiple  sub-pixel  locations  of  objects  are  used  to  create  modeled  PSFs  to  compare 
potential  objects  against.  A  key  factor  in  the  development  of  the  MHT  is  assigning  the  costs 
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and  prior  probabilities  of  the  hypotheses,  which  represent  two  of  the  research  questions 
and  are  investigated  in  Chapters  3  and  4  respectively.  Figure  2.4  shows  the  structure 
of  the  relevant  detection  algorithms  described  in  this  section.  As  shown  in  Figure  2.4, 
the  detection  algorithms  considered  to  this  point  can  be  described  in  a  hierarchy.  This 
dissertation  focuses  on  the  matched  filter  MHT,  and  the  choices  for  assigning  costs  and 
prior  probabilities. 


Figure  2.4:  Hierarchy  of  space  object  detection  techniques  considered  in  this  research. 
This  visualization  demonstrates  how  the  point  detection,  matched  filter,  BHT,  and  MHT 
algorithms  are  related. 


In  addition,  the  choice  for  the  received  noise  can  be  applied  to  any  of  the  detection 
algorithms  shown  in  Figure  2.4,  potentially  changing  the  performance  of  the  algorithms. 
One  key  factor  in  building  a  detector  is  accurately  modeling  the  distribution  of  the  expected 
received  data.  In  Chapters  3  and  4,  it  is  assumed  that  the  received  data  is  Gaussian,  but  a 
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method  is  presented  for  implementing  non-Gaussian  received  data  tests  with  different  costs 
and  prior  probability  assumptions.  In  Chapter  5,  a  new  model  for  noise  in  the  PSF  is 
developed,  and  a  detection  algorithm  based  on  the  model  is  used. 

The  MHT  proposed  in  this  research  is  based  on  a  Bayes  risk  described  by  Van  Trees 
[25].  In  Bayes  risk  there  are  two  important  variables  that  can  be  assigned  that  change 
the  form  and  effectiveness  of  the  test.  They  are  the  costs  associated  with  selecting  each 
hypothesis  and  the  prior  probabilities,  or  the  probability  of  each  hypothesis  occurring. 
While  multiple  hypothesis  detection  based  on  Bayes  risk  is  a  relatively  new  research  area 
in  SDA,  it  has  applications  in  several  other  areas.  One  of  these  is  the  detection  and 
classification  of  reflecting  objects  in  radar  imagery  investigated  by  Ertin  [26].  The  authors 
present  an  unequal  cost  MHT,  by  penalizing  misclassification  of  objects  differently  than 
missed  detections  or  false  alarms.  A  key  difference  is  that  the  choice  for  the  cost  of  a 
misclassification  still  requires  the  algorithm  to  make  a  classification  decision,  which  is  not 
done  in  Chapter  3.  Another  example  of  unequal  cost  MHT  is  the  removal  of  symbols 
with  errors  in  a  frequency-hop  communication  system  investigated  by  Baum  and  Pursley 
[27].  Again,  different  choices  for  costs  and  received  data  assumptions  lead  to  different 
algorithms. 

Unequal  prior  probability  has  use  in  other  research  areas  outside  of  SDA  and  space 
object  detection.  Often  in  cases  where  the  input  conditions  might  change,  and  adapting 
the  assumptions  or  inputs  to  the  algorithm  may  provide  additional  performance  gains. 
These  input  conditions  are  the  prior  probabilities.  Cognitive  radar  [28],  neural  networks 
[29]  and  adaptive  algorithms  [30,  31],  Bayes  estimators  [32],  and  quantization  of  prior 
probabilities  [33]  are  research  areas  that  have  investigated  this  effect.  In  some  cases  the 
focus  is  comparing  equal  vs.  unequal  prior  probabilities  only,  but  some  research  has  also 
been  done  on  a  non-constant  assignment  of  priors.  In  the  case  of  the  SST  and  SDA,  the 
priors  may  change  between  collections  or  frames  due  to  changing  atmospheric  conditions. 
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Decisions  about  how  often  the  priors  need  to  be  updated,  or  the  optimal  feedback  for  the 
priors  are  not  covered  in  this  research.  In  the  next  section,  the  relevant  background  for  the 
optical  model  is  discussed. 

2.2  Optical  Modeling  Overview 

A  key  element  in  this  research  is  choosing  an  accurate  optical  model.  The  optical 
model  developed  here  is  used  for  several  purposes  throughout  this  research,  including 
simulating  data  and  creating  the  PSFs  for  the  matched  filter  algorithms.  This  section 
provides  a  brief  overview  of  the  model  and  details  where  the  model  is  described  further. 

Depending  on  the  research  area,  different  aspects  of  the  optical  model  are  used.  The 
main  components  of  an  optical  model  can  be  described  by  the  following  four  factors:  the 
propagation  model,  the  model  for  aberrations,  atmospheric  model,  and  the  received  noise 
model.  In  this  dissertation,  the  light  captured  by  the  telescopes  is  generally  incoherent. 
This  is  because  the  source  of  the  light  from  RSOs  in  orbit  is  reflected  incoherent  sunlight. 
Additional  complications  and  considerations  are  required  for  coherent  light,  but  these 
topics  are  not  considered  in  this  research. 

The  first  step  in  creating  an  accurate  optical  model  is  recreating  the  electromagnetic 
field  as  it  propagates  from  one  point  in  space  to  another.  Rayleigh-Sommerfeld  propagation 
is  derived  from  Maxwell’s  equation  and  describes  an  electromagnetic  field  as  the  result 
of  a  wave  traveling  from  one  plane  to  another.  Although  it  is  the  most  accurate  way  to 
completely  describe  a  field  after  a  propagation,  it  is  also  the  most  computationally  complex 
method  for  propagating  a  field.  For  the  large  CCD  arrays  considered  in  this  research,  it  is 
not  feasible.  Rayleigh-Sommerfeld  propagation  can  be  reduced  by  making  an  assumption 
about  the  distance  from  the  source  to  the  receiver  plane,  along  with  the  relative  size  of 
source  and  receiver  planes  [34].  These  factors  determine  which  method  can  be  accurately 
implemented. 
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(2.1) 


z3>>£f[(x“^)2+(2'“’')2]l 

z  is  propagation  distance  between  the  source  and  received  plane,  A  is  the  wavelength 
being  considered,  x,  y  are  receiver  plane  coordinates  and  rj  are  source  plane  coordinates. 
If  equation  (2.1)  is  true,  a  Fresnel  propagation  can  be  implemented.  A  Fresnel  propagation 
is  less  computationally  complex  than  the  Rayleigh-Sommerfeld  propagation.  The  Fresnel 
approximation  can  be  further  reduced  to  a  Fraunhofer  propagation.  Fraunhofer  can  be 
thought  of  as  the  far  field  effect  of  a  propagation,  or  in  the  case  where  a  focusing  optic  is 
employed.  It  is  accurate  when  the  following  condition  is  met: 
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where  k  is  the  wave  number  =f.  If  this  condition  is  met,  it  can  be  considered  valid  to 
utilize  a  Fraunhofer  propagation.  The  benefit  of  using  a  Fraunhofer  propagation  is  that  it 
can  be  implemented  with  a  two-dimensional  Fourier  Transform. 


U  (x,  y) 


gjkZgj^i^+y2) 


jAz 


oo 


rj)exp 


2n 

-J-rix^  +  yri) 

Az 


d^di] 


(2.3) 


C/(£,  rj)  is  the  field  at  the  source  plane,  and  U(x,y )  is  the  resulting  field  at  the  receiver 
plane.  Throughout  this  dissertation,  Fraunhofer  propagation  is  used.  To  generate  the  PSF, 
the  source  field  is  assumed  to  be  a  single  point  source  located  at  a  significant  distance  from 
the  aperture.  The  result  of  this  first  long  propagation  is  a  plane  wave  normal  to  the  optic 
axis  of  the  telescope  system.  The  second  propagation  needed  is  from  the  telescope  aperture 
to  the  CCD  received  plane.  This  is  where  the  Fraunhofer  propagation  is  computed.  The 
propagation  from  aperture  to  receiver  needs  to  take  into  account  other  factors,  including 
the  telescope  aberrations  and  atmosphere,  which  are  discussed  next. 

An  imaging  system  can  be  considered  diffraction  limited  “if  a  diverging  spherical 
wave,  emanating  from  a  point-source  object,  is  converted  by  the  system  into  a  new  wave, 
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again  perfectly  spherical,  that  converges  toward  an  ideal  point  in  the  image  plane  [34].”  In 
other  words,  a  diffraction  limited  system  is  an  ideal  system  that  contains  no  imperfection 
on  the  lenses  or  mirrors.  The  diffraction  limited  system  has  the  largest  spatial  frequency 
content  allowed  by  the  aperture  dimensions.  Any  aberrations  added  to  a  diffraction  limited 
system  limit  the  maximum  allowable  frequency. 

In  a  non-diffraction  limited  system,  the  defects  in  the  optics  used,  or  aberrations,  affect 
the  formation  of  the  image  and  generally  need  to  be  accounted  for  to  create  an  accurate 
model.  In  Chapters  3  and  4,  these  aberrations  are  included  in  the  model  and  described. 
In  Chapter  5  aberrations  are  not  included.  This  is  because  the  focus  of  the  research  is  on 
the  temporal  changes  in  the  PSF  caused  by  the  atmosphere,  and  the  lens  aberrations  are 
assumed  to  not  change  temporally  over  the  integration  period  of  the  system. 

The  effect  of  the  atmosphere  on  the  optical  model  is  captured  in  two  different  methods 
in  this  dissertation.  These  methods  are  dependent  on  the  length  of  integration  time  used  in 
the  data  collection.  The  two  methods  are  a  long  exposure  atmosphere  and  a  short  exposure 
atmosphere  model.  More  details  on  the  model  for  the  atmosphere  are  discussed  in  section 
2.4. 


2.3  Sampling  Theory 

The  primary  motivation  for  using  a  MHT  is  to  overcome  the  effects  of  undersampling 
by  the  CCD  pixels.  If  a  telescope  system  is  spatially  undersampled,  a  small  shift  of  where 
an  object  is  formed  within  a  single  pixel  can  result  in  a  different  shape  or  distribution  of 
intensity.  Instead  of  shifting  and  retaining  the  spatial  information,  an  aliased  PSF  loses 
spatial  information  and  can  have  a  different  shape  completely.  This  presents  a  problem  in 
matched  filter-based  detection  algorithms,  where  the  goal  is  to  match  objects  in  the  data 
to  the  hypothesized  PSFs.  Figure  2.5  shows  four  examples  of  an  undersamped  PSF  model 
generated  with  a  small  sub-pixel  shifts. 
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Figure  2.5:  Example  of  four  simulated  PSFs  impacted  by  spatial  undersampling  caused  by 
CCD  pixels.  Each  PSF  is  shifted,  but  also  changed  in  the  shape  and  intensity  distribution. 


As  shown  in  Figure  2.5,  the  intensity  and  distribution  of  the  PSF  changes  as  the  RSO 
is  moved  within  a  pixel.  If  the  RSO  is  properly  sampled,  the  PSFs  would  have  the  same 
shape,  and  only  be  shifted  slightly  in  position.  A  common  method  of  determining  the 
sampling  required  is  using  angular  resolution  found  through  a  Rayleigh  criterion.  In  [23], 
O’Dell  shows  that  in  a  matched  filter  detection  algorithm,  Rayleigh  always  undersamples 
when  compared  to  the  required  Nyquist  sampling.  The  required  Nyquist  sampling  is  found 
by  combining  the  maximum  frequency  present  in  a  diffraction  limited  system  with  the 
Nyquist  theorem  on  minimum  required  sampling.  This  gives  the  following  required  spatial 
sampling,  A.s,  described  by  Goodman  [34]: 


As  <  —  (2.4) 

where  A  is  the  center  wavelength  observed,  /  is  focal  length  of  the  system,  and  d  is  the 
diameter  of  the  pupil.  Fooking  at  the  variables  of  equation  (2.4),  it  is  evident  that  there  is 
a  limited  amount  of  user  control  in  the  required  Nyquist  sampling.  The  center  wavelength 
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is  a  property  of  the  observed  light.  The  main  method  of  selecting  the  desired  sampling 
rate  is  the  ratio  of  4,  but  these  parameters  are  essential  to  other  portions  of  system  design. 
These  design  components  include  the  Field  of  View  (FoV)  of  the  telescope  and  the  physical 
dimensions  of  the  system. 

Substituting  the  SST  parameters  into  equation  (2.4),  a  center  wavelength  of  500nm 
with  a  focal  length  and  pupil  diameter  of  3.5m,  gives  a  required  sampling  size  of  0.25/mi. 
The  SST  pixels  are  15/rni  square,  but  are  grouped  2x2  in  these  data  collections.  The 
binned  pixels  give  an  actual  sampling  size  of  30/iin.  The  difference  between  required 
and  actual  sampling  in  this  case  is  120  times.  If  the  system  is  not  diffraction  limited,  the 
actual  undersampling  factor  will  be  much  less.  This  is  due  to  the  pupil  diameter  d  being 
limited  by  the  effective  seeing  parameter,  ra.  Seeing  parameter  values  change  depending 
on  atmospheric  conditions,  but  are  typically  smaller  than  10cm  at  the  Socorro,  NM  site. 

2.4  Developing  an  Accurate  Model  for  Scintillation  Noise  Present  in  a  PSF 

The  noise  present  in  an  image  captured  through  a  telescope  can  have  several  potential 
sources  including  read  out  noise,  dark  current,  photon  counting  noise,  and  noise  caused 
by  the  atmosphere.  Traditionally,  SDA  detection  algorithms  have  utilized  a  Central  Limit 
Theorem  justification  to  assume  that  the  combination  of  all  of  the  noise  sources  present 
in  the  system  result  in  an  overall  Gaussian  distribution  of  noise  within  the  received  data. 
Multiple  research  efforts  have  looked  at  detection  algorithms  based  on  the  noise  in  the 
received  data  not  being  Gaussian.  In  multiple  research  efforts  it  is  assumed  that  the 
noise  is  dominated  by  Poisson  noise,  and  other  sources  do  not  contribute  significantly 
[22,  35].  In  these  cases,  a  more  complicated  relationship  is  found  to  be  necessary  in  order  to 
implement  the  algorithm.  One  example  of  building  a  detector  optimized  for  Poisson  noise 
is  investigated  by  Pohlig  [22].  In  this  case,  the  dominant  source  of  noise  in  the  received  data 
is  the  random  arrival  time  of  photons.  The  author  states  this  can  be  achieved  through  a  very 
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low  noise  CCD  in  low  intensity  collection.  In  [35],  Peterson  assumes  that  the  distribution 
of  the  received  data  is  Poisson,  and  the  author  achieves  similar  results. 

The  goal  in  this  portion  of  research  is  to  account  for  the  two  phenomenon  present, 
and  combine  them  into  a  single  model.  These  phenomenon  are  the  atmosphere  and  photon 
counting  noise.  It  is  well  understood  that  noise  due  to  random  photon  arrival  times  follows 
a  Poisson  distribution  [17].  The  impact  of  the  atmosphere,  and  the  fluctuations  it  causes 
in  received  intensity  is  less  understood.  This  section  outlines  atmospheric  theory  that  has 
been  developed,  as  well  as  current  methods  for  overcoming  its  effects. 

2.4.1  Overview  of  Atmosphere  Theory. 

The  atmosphere  causes  random  wavefront  errors  to  any  field  that  propagates  through 
it,  and  these  wavefront  errors  are  due  to  changes  in  temperature  and  wind  at  different 
points  in  the  propagation  path  that  affect  the  index  of  refraction  [36].  The  degradation 
in  image  quality  due  to  the  atmosphere  presents  significant  challenges  across  several  types 
of  imaging  systems.  Both  coherent  and  incoherent  optical  systems,  from  laser  systems 
and  communication  to  ground  based  telescope  imaging  account  for  the  effects  of  the 
atmosphere. 

These  systems  utilize  different  techniques  to  improve  the  performance  of  the  optical 
system.  Given  the  fact  that  the  atmosphere  behaves  as  a  random  process,  it  is  impossible  to 
completely  describe  its  behavior  with  a  deterministic  set  of  equations.  Instead,  many  efforts 
have  looked  at  determining  the  statistics  of  the  atmosphere.  Assuming  the  atmosphere 
is  statistically  homogeneous  and  isotropic,  the  covariance  can  be  used  to  give  a  Power 
Spectral  Density  (PSD)  as  detailed  by  Andrews  and  Phillips  [36].  The  PSD  is  a  function 
of  the  strength  of  turbulence  in  the  atmosphere,  C2N.  In  horizontal  propagation,  C2N,  is  often 
assumed  to  be  constant.  In  vertical  or  slant  imaging,  as  in  stellar  observations,  C2N  is  often 
modeled  as  a  function  of  elevation,  h.  Figure  2.6  shows  an  example  plot  of  a  popular  model 
for  C2n,  the  Hufnagel- Valley. 
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Figure  2.6:  Hufnagel- Valley  C2W  as  a  function  of  elevation  in  meters.  C\,  represent  the 
strength  of  the  turbulence  in  the  atmosphere. 


The  Hufnagel- Valley  profile  depends  on  the  value  of  turbulence  on  the  ground,  A. 
This  value  affects  the  profile  up  to  approximately  3000m.  The  value  A  =  1.7xl0-14  is 
selected  for  the  commonly  used  H-Vs  model.  Combining  the  statistical  description  of  the 
atmosphere  with  known  wave  propagation  equations,  theoretical  derivations  can  be  made 
for  the  received  fields,  intensity,  and  intensity  fluctuations.  These  are  related  to  the  first, 
second,  and  fourth  order  moments  respectively. 

The  random  phase  errors  in  light  propagated  through  the  atmosphere  cause  changes 
in  the  shape,  brightness,  and  location  of  the  image  formed  at  the  focal  length  or  CCD 
of  an  imaging  system.  These  changes  manifest  as  intensity  fluctuations  in  the  received 
image.  These  fluctuations  are  known  as  scintillation.  Scintillation  can  present  as  a 
temporal  fluctuation,  such  as  a  star  twinkling  over  time,  or  as  a  spatial  fluctuation,  known 
as  speckle  [36].  Characterizing  how  scintillation  can  impact  the  PSF  will  result  in  an 


21 


improved  statistical  description  of  noise  in  the  received  data.  Knowledge  of  noise  statistics 
could  potentially  contribute  to  improving  detection  algorithm  performance.  A  thorough 
summary  of  the  progress  towards  understanding  and  correcting  for  the  earth’s  atmosphere 
is  described  by  Roddier  in  [37]. 

Accounting  for  the  effect  of  the  atmosphere  on  the  PSF  is  dependent  on  the  exposure 
time,  Ts,  used.  In  the  short  exposure  atmosphere,  the  intensity  has  less  spread,  but  contains 
more  variance  in  shape,  intensity,  and  location.  The  long  exposure  atmosphere  averages 
the  random  phase  effects,  and  as  a  result  deceases  the  fluctuations,  causing  a  larger  spatial 
spread  in  the  PSF.  For  exposure  times  much  greater  than  10ms,  atmospheric  effects  are 
effectively  averaged,  and  a  long  exposure  atmosphere  model  can  be  accurately  used.  In  the 
long  exposure  atmosphere  model,  the  received  data  is  assumed  to  be  constant  in  time,  and 
is  a  function  of  seeing,  wavelength,  focal  length,  and  position.  This  is  due  to  the  fact  that 
the  intensity  variance  over  time,  aj,  effectively  goes  to  zero.  For  Ts  much  less  than  10ms, 
the  effects  of  the  atmosphere  are  “frozen  [17].”  This  leads  to  a  large  variance,  07.  Figure 
2.7  demonstrates  two  simulated  images  captured  with  a  short  exposure  and  long  exposure 
time. 


Figure  2.7:  Simulated  star  captured  with  both  (A)  short  and  (B)  long  exposure  time,  Ts. 
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In  (A),  the  short  exposure  image  has  two  factors  impacting  its  appearance.  The  first 
is  that  since  a  small  exposure  time  is  used,  the  captured  intensity  is  quite  low.  As  a  result, 
it  is  more  difficult  to  separate  the  object  from  background  noise.  In  addition,  the  object 
may  contain  spatial  variations  due  to  the  atmosphere.  If  a  second  short  exposure  image 
is  shown,  the  effect  of  temporal  scintillation  could  also  be  demonstrated  by  the  change  in 
intensity  between  the  two  images.  In  (B),  the  long  integration  time  has  allowed  for  more 
intensity,  distinguishing  the  object  more  clearly  from  the  background.  In  addition,  the 
spatial  dimensions  of  the  object  have  increased  due  to  the  averaging  of  the  atmosphere. 

An  example  of  a  system  operating  in  the  short  exposure  model  is  found  in  the  laser 
communication  area.  Characterizing  intensity  fluctuations  are  essential  to  the  field  of  laser 
communications.  The  atmospheric  degradation  effect  on  laser  beams  causes  a  loss  of 
intensity  and  the  spreading  of  the  laser  spot.  Due  to  the  high  data  rates  typically  used, 
these  systems  operate  in  the  short  exposure  regime.  A  key  factor  in  a  laser  communication 
system  is  signal  reliability,  which  is  related  to  the  PDF  of  the  intensity  according  to 
Al-Habash  et.  al.  [38].  For  this  reason,  several  research  efforts  have  investigated  the 
fundamental  distribution  of  noise  in  the  propagation  of  light  through  the  atmosphere. 
Alternatively,  a  telescope  imaging  distant  objects  would  need  a  large  time  period  to  collect 
a  sufficient  number  of  photons  from  the  dim  object.  Identifying  the  characteristics  of 
different  atmospheric  effects  is  important,  and  can  help  to  determine  appropriate  methods 
for  correcting  them.  Removing  atmospheric  spread  provides  an  estimate  of  the  original 
object,  increasing  clarity  and  resolution.  More  information  about  these  efforts  is  included 
in  the  next  section. 

2.4.2  Methods  for  correcting  atmospheric  effects. 

There  are  several  ways  to  mitigate  atmospheric  effects.  Adaptive  Optics  (AO)  and 
image  restoration  are  two  popular  methods.  They  can  be  used  individually  or  together  to 
achieve  increased  performance.  In  AO,  the  wavefront  errors  are  estimated,  and  a  control 
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system  is  used  to  actuate  deformable  mirrors  to  correct  for  these  errors  before  the  image  is 
formed  at  the  CCD.  Adaptive  optics  are  implemented  during  the  design  phase  of  system 
development.  Using  AO,  a  portion  of  the  wavefront  errors  are  corrected  before  the  image 
is  digitally  captured.  This  provides  resolution  benefits  without  post-processing  of  the  data. 
The  additional  performance  comes  at  a  cost  in  hardware,  sensing  and  control,  and  cost. 
Depending  on  the  desired  application  and  the  benefits  provided,  these  additional  costs  are 
not  always  justifiable.  In  AO  corrected  systems,  there  are  often  residual  phase  errors  that 
are  not  corrected  for.  To  correct  for  remaining  errors,  additional  methods  such  as  image 
restoration  can  still  be  used.  In  [39,  40],  the  issue  of  image  processing  in  AO  corrected 
images  is  investigated  by  Racine  et.  al.  and  Fusco  et.  al. 

Adaptive  optics  systems  are  not  investigated  further  in  this  research  for  two  primary 
reasons.  One  of  these  reasons  is  that  the  SST  system,  which  is  the  focus  of  this  research, 
does  not  utilize  an  AO  system.  AO  is  not  utilized  in  SST  mainly  due  to  the  large  FoV 
[41].  AO  operates  on  the  assumption  that  all  portions  of  the  received  wave  propagate 
through  the  same  atmosphere  or  one  that  is  highly  correlated.  In  a  large  FoV  collection, 
this  assumption  is  not  valid.  Current  research  efforts  have  demonstrated  AO  correction 
over  a  2  arcmin  FoV —  [42],  a  small  percentage  of  the  FoV  of  SST.  The  second  reason  AO 
are  not  investigated  further  is  that  the  techniques  developed  here  can  still  be  utilized  in  an 
adaptive  optics  system.  Figure  2.8  shows  an  example  block  diagram  of  an  adaptive  optics 
telescope  system. 

A  second  method  often  utilized  is  correcting  for  the  atmospheric  effects  with  image 
processing.  In  this  case,  no  additional  hardware  is  needed.  The  image  formed  at  the 
CCD  still  contains  the  effects  of  the  atmosphere,  and  post-processing  is  used  to  obtain 
the  original  image  or  to  account  for  the  atmosphere.  Within  the  field  of  image  processing, 
there  are  two  important  cases  that  lead  to  different  approaches.  These  cases  are  resolved 
and  unresolved  imaging.  In  a  resolved  image,  separate  portions  of  an  imaged  object  are 
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Figure  2.8:  Block  diagram  of  an  AO  telescope  system  [43]. 


able  to  be  differentiated  from  each  other.  In  an  unresolved  image,  the  object  may  occupy  a 
single  or  small  number  of  pixels,  and  much  if  not  all  spatial  information  about  the  object 
is  lost.  The  ability  to  resolve  an  object  is  a  function  of  several  factors,  including  CCD 
sampling,  lenses,  wavelength,  and  the  distance  and  size  of  the  target.  In  resolved  imaging, 
it  is  often  desirable  to  use  signal  processing  to  remove  the  effects  of  the  atmosphere.  These 
systems  are  typically  trying  to  capture  physical  characteristics  of  the  object  being  imaged, 
and  the  blurring  caused  by  the  atmosphere  degrades  and  limits  this  capability. 

There  are  several  proposed  techniques  for  image  restoration  in  resolved  imaging.  One 
common  practice  is  to  use  deconvolution  [44,  45].  In  deconvolution,  an  estimate  of  the 
truth  image  is  computed  through  the  reversal  of  the  process  of  image  formation.  From 
optical  theory,  an  image  is  formed  through  the  convolution  of  the  truth  object  with  the  PSF 
of  the  optical  system  [34].  Using  knowledge  about  the  possible  and  likely  PSF,  an  estimate 
can  be  formed  of  the  truth  object.  Additionally,  there  are  methods  for  estimating  both  the 
PSF  and  the  original  object,  known  as  blind  deconvolution. 

In  an  unresolved  imaging  system,  the  spatial  information  of  the  object  may  not  be  as 
important.  As  a  result,  it  is  not  necessary  to  implement  costly  image  restoration  algorithms. 
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One  example  of  where  this  type  of  imaging  may  occur  is  in  space  imaging  of  distant  or 
small  objects,  as  is  the  case  in  this  research.  Although  there  may  be  instances  of  image 
restoration  in  space  imaging,  including  attempting  to  resolve  binary  star  systems,  this  effort 
does  not  use  these  techniques. 

In  an  optical  system  built  for  detection  of  space  objects,  the  algorithm  used  to  detect 
the  objects  is  key.  In  this  situation,  the  objective  is  to  create  a  binary  map  of  object 
locations  from  an  image.  A  detection  algorithm  is  used  to  determine  if  each  pixel  is  one 
or  zero  depending  on  if  an  object  is  considered  to  be  in  that  pixel  or  not.  This  map  of  the 
detected  object  is  then  used  for  further  image  processing  techniques  to  determine  additional 
information  for  cataloging  objects.  Such  information  may  include  orbit  determination, 
size,  and  comparison  against  already  known  objects.  As  described  earlier,  these  algorithms 
depend  on  knowledge  of  the  distribution  of  the  noise  in  the  image.  In  these  situations, 
the  noise  caused  by  the  atmosphere  can  be  accounted  for  by  understanding  the  statistical 
fluctuations  it  causes. 

This  type  of  direct  detection  system  can  be  useful  in  a  situation  where  large  amounts 
of  data  need  to  be  processed.  The  image  restoration  methods  described  previously  can 
be  computationally  expensive.  In  a  large  survey  telescope  like  the  SST,  the  required 
processing  times  are  not  feasible.  In  addition,  the  spreading  caused  by  the  atmosphere 
helps  create  distinct  PSF  shapes  that  can  be  used  for  match  filter  detection  algorithms, 
which  are  discussed  in  Chapter  3.  This  research  focuses  on  the  analysis  of  the  distribution 
of  the  received  data  to  directly  perform  detection  algorithms  without  AO  or  the  image 
restoration  techniques  described  above.  There  are  several  approaches  for  obtaining  a  noise 
distribution.  These  methods  are  described  in  the  next  section. 

2.4.3  Statistical  models  for  scintillation. 

There  are  three  practical  methods  for  developing  a  statistical  model  for  the  noise  in 
the  PSF.  These  methods  are  a  theoretical  based  approach,  simulation,  and  analysis  of 


26 


measured  data.  Each  of  these  methods  present  benefits  and  challenges.  If  the  results  derived 
from  these  methods  agree,  it  increases  confidence  in  the  achieved  result  from  any  of  the 
individual  methods. 

Determining  a  statistical  model  from  theory  is  the  most  robust  method  for  identifying 
the  true  distribution.  This  approach  is  independent  of  the  telescope  or  imaging  system 
used,  and  is  therefore  applicable  in  any  optical  system.  One  drawback  is  that  due  to 
its  complicated  nature,  it  may  not  be  realistic  to  develop  a  closed  form  solution  to  this 
problem.  Starting  from  the  first  principles  of  atmospheric  theory  is  the  most  accurate 
way  to  build  an  improved  detector.  Several  research  efforts  have  attempted  to  determine 
distributions  of  received  intensity  through  the  atmosphere.  These  efforts  have  focused  on 
both  characterizing  the  statistics  of  the  atmosphere  [46-48]  and  the  resulting  PSF  formed 
from  an  optical  system  viewing  through  the  atmosphere  [49,  50].  One  generally  accepted 
model  of  the  intensity  fluctuations  due  to  the  atmosphere  is  log-normal  distribution.  The 
log-normal  is  derived  from  theory  developed  by  Tatarskii  in  [51,  52],  and  is  experimentally 
verified  in  multiple  physical  experiments  [53,  54].  The  drawback  to  these  methods  is  that 
they  assume  that  the  noise  caused  by  the  atmosphere  is  the  only  significant  source,  and 
therefore  is  the  dominant  source  of  noise  in  the  image.  This  may  be  the  case  for  certain 
intensity  objects,  atmospheric  conditions,  or  exposure  times.  However,  this  assumption  is 
not  made  in  this  research. 

A  second  method  for  determining  a  statistical  model  is  to  collect  and  analyze  a  large 
amount  of  data  to  develop  a  model.  This  approach  does  not  depend  on  deriving  the 
fundamental  properties  of  the  atmosphere,  but  attempts  to  characterize  it  through  its  impact 
on  measured  data.  The  drawback  of  this  method  is  that  the  results  are  dependent  on  a  large 
number  of  data  collections,  and  the  analysis  of  those  collections.  There  are  several  relevant 
examples  of  published  research  using  this  method.  This  observational  method  is  used  by 
Jee  et.  al.  in  [55],  but  the  authors  characterize  spatial  fluctuations  in  the  PSF  rather  than 
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temporal.  Additionally,  this  work  focuses  on  the  Advanced  Camera  for  Surveys  (ACS) 
which  is  installed  on  the  Hubble  Space  Telescope  (HST).  This  is  a  key  difference  between 
this  research  effort  and  [55].  Since  HST  is  outside  of  the  Earth’s  atmosphere,  the  PSF 
fluctuations  are  not  caused  by  the  atmosphere,  but  by  other  factors. 

Another  research  effort  that  investigates  the  statistical  and  temporal  properties  of 
scintillation  is  completed  by  Dravins  et.  al.  in  [56].  In  this  paper,  two  statistical  aspects  of 
intensity  fluctuations  are  investigated.  The  first  uses  a  photon  counter  to  determine  a  density 
function  for  the  number  of  received  photons.  The  second  statistical  aspect  investigated  is 
the  autocorrelation.  The  autocorrelation  provides  insight  into  the  behavior  of  the  intensity 
fluctuations  over  time.  Due  to  the  nature  of  the  data  collection  method  described  in  this 
paper,  the  authors  investigate  short  exposure  times,  in  the  sub-millisecond  regime.  It  is 
shown  that  over  a  specified  time  separation,  the  correlation  drops  to  a  negligible  amount. 
At  this  time  separation,  it  can  be  assumed  that  the  two  images  are  independent.  This  fact 
could  be  used  to  characterize  longer  exposure  times  that  might  be  utilized  by  the  SST.  To 
match  the  collected  data  to  potential  distribution  fits,  two  methods  may  be  used.  The  first 
method  is  to  fit  histograms  of  the  received  data  to  PDF  and  utilize  an  error  metric  like  Root 
Mean  Squared  Error  (RMSE).  The  second  method  is  to  investigate  the  statistical  moments 
of  the  received  data.  With  this  method,  the  authors  experimentally  determine  that  a  beta 
distribution  of  the  second  kind  provides  the  closest  match  to  the  statistical  moments  under 
their  assumptions.  The  authors  have  no  physical  motivation  for  choosing  this  distribution, 
only  stating  that  it  provides  the  closest  match  to  the  third  and  fourth  order  moments  of  the 
distributions  tested. 

A  third  approach  to  characterize  the  noise  present  is  to  utilize  a  simulation.  The  benefit 
of  a  simulation  is  that  a  large  amount  of  data  can  be  easily  generated  for  little  to  no  cost. 
This  can  solve  some  of  the  problems  presented  by  the  measured  data  method.  A  drawback 
to  this  method  is  that  several  limiting  assumptions  are  made  to  make  the  simulation  feasible. 
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In  addition,  both  measured  data  and  simulation  depend  on  utilizing  specific  telescope 
parameters  for  several  properties,  including  lens  aberrations  and  telescope  diameter  to 
atmospheric  seeing  ratio.  Performing  a  simulation  relies  on  modeling  two  key  components 
of  the  system:  the  atmosphere  and  the  optical  system.  Building  an  optical  model  is 
straightforward  and  well  understood.  Only  a  few  assumptions  such  as  a  Fraunhofer 
propagation  can  be  used.  Simulating  light  propagating  through  the  atmosphere  presents 
several  challenges.  The  first  is  whether  to  represent  the  atmosphere  as  a  single  or  multiple 
phase  screen.  Depending  on  assumptions  of  the  viewing  angle,  either  can  be  justifiable. 

In  this  research,  a  single  random  phase  screen  is  used,  and  its  details  are  developed 
in  Chapter  5.  In  addition,  there  are  several  proposed  methods  for  creating  random  phase 
screens.  They  are  a  Fourier  Transform-based  phase  screen  described  by  Schmidt  [57]  and 
a  Zemike  polynomial  based  phase  screen  developed  by  Cain  and  Richmond  [58].  Both  of 
these  methods  use  the  underlying  statistics  developed  in  Kolmogorov  theory  [59,  60],  but 
implement  these  statistics  in  different  ways.  In  [61],  Roddier  et.  al.  show  that  the  wavefront 
errors  due  to  the  atmosphere  can  be  accurately  represented  with  Zemike  polynomials. 
Using  the  models  described  above,  the  simulated  intensity  can  be  observed  as  a  function  of 
both  exposure  time  and  seeing  conditions.  From  this  data,  a  statistical  fit  can  be  made  to 
determine  which  distribution  most  closely  matches  the  simulated  results. 

This  dissertation  will  first  theoretically  derives  the  solution  for  the  distribution  of 
intensity  fluctuations  in  a  received  PSF.  To  provide  a  validation  of  the  theory,  simulated 
space  object  data  will  be  analyzed.  They  will  provide  solutions  in  the  case  where  an  analytic 
solution  is  not  feasible.  Chapter  3  will  outline  the  development  of  an  improved  Multiple 
Hypothesis  detection  algorithm  based  on  a  Gaussian  assumption  of  noise  in  the  collected 
image. 
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III.  Improved  Multiple  Hypothesis  Test  through  an  Unequal  Cost  Assumption 


In  this  Chapter  a  new  detection  algorithm  that  could  be  utilized  in  SDA  telescopes 
is  proposed.  Section  3.1  presents  an  optical  model  for  the  SST  system  that  is  used  in  the 
algorithm.  Section  3.2  describes  the  theory  involved  in  developing  the  detection  algorithm. 
Section  3.3  discusses  the  experimental  setup  used  to  collect  data  for  this  Chapter.  Section 
3.4  outlines  results  and  data  analysis.  Section  3.5  discusses  a  method  for  full  frame 
implementation  of  the  detection  algorithm.  Section  3.6  discusses  the  conclusions  made 
based  on  the  presented  results  and  potential  future  work. 

3.1  SST  Optical  Model 

The  image  created  by  an  optical  system  viewing  a  point  source  or  spatial  impulse 
is  known  as  the  PSF.  The  algorithm  developed  for  this  chapter  requires  a  model  for  the 
expected  PSF  of  the  objects  being  investigated.  The  SST  views  space  objects  that  are  either 
relatively  small  and  located  in  earth  orbit,  or  very  distant  objects  like  stars  and  asteroids. 
Using  geometric  optics,  the  ratio  of  pixel  size  to  focal  length  can  be  compared  to  the  size 
of  the  object  and  distance  from  the  SST.  Using  this  relation,  it  is  evident  that  all  of  these 
potential  objects  are  effectively  point  sources  to  the  telescope. 

As  described  in  Chapter  2,  if  the  conditions  are  met,  a  Fourier  Transform  can  be  used 
to  perform  a  field  propagation.  Goodman  demonstrated  that  the  PSF  of  an  optical  system, 
h0pt(x,y),  can  be  found  by  the  following  relation  [34]: 

hopt(x,y)  =  \T{P(m,n)}\2  (3.1) 

where  x,y  are  spatial  distance  pixel  coordinates  in  the  detector  plane,  and  m,n  are  pixel 
coordinates  in  the  pupil  plane.  P(m,  n )  is  a  pupil  function  that  mathematically  describes  the 
effect  of  the  pupil  on  incoming  light,  and  T  is  a  two-dimensional  Fourier  transform.  This 
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relation  holds  for  a  flat  focal  plane  array,  which  the  SST  does  not  have.  It  is  assumed  that 
the  SST  does  not  have  a  focal  plane  array  in  this  research  effort.  This  implies  there  may  be 
small  errors  in  PSFs  created  off  the  optic  axis.  The  creation  of  optimal  PSF  is  not  the  focus 
of  the  research,  and  all  algorithms  tested  use  the  same  flat  focal  plane  array  assumption. 

In  a  diffraction  limited  optical  system  consisting  of  a  perfect  lens  or  mirror,  the  pupil 
function  consists  of  only  the  geometry  of  the  pupil,  P(m,  n)  =  A(m ,  n),  and  contains  no 
other  phase  distortions.  A(m,  n )  is  an  amplitude  function  that  is  one  or  zero,  depending  on 
if  light  is  able  to  pass  through  the  pupil  at  the  m,  n  pixel  location.  In  this  optical  model, 
the  SST  is  assumed  to  have  a  3.5m  primary  mirror  and  a  1.75m  obscuration.  The  physical 
telescope  has  secondary  mirror  arms  and  other  minor  obstructions,  but  these  objects  do  not 
have  a  significant  impact  on  the  produced  image. 

In  a  non-ideal  imaging  system,  imperfections  in  the  lenses  or  mirrors  cause  phase 
distortions  to  any  light  passing  through  the  optics.  These  distortions  are  modeled  as  phase 
fluctuations  to  the  pupil  function  [34]: 

P(m,  n )  =  A(nu  n)  exp  [j0o(m,  ri)\  (3.2) 

The  phase  aberrations  90(m,n),  are  expressed  as  the  sum  of  a  set  of  orthonormal 
Zernike  polynomials,  with  each  polynomial  representing  a  type  of  phase  distortion  [62]. 

9 aim,  n)  =  ^  akZk{m ,  n)  (3.3) 

k 

The  weighting  coefficients  ak  represent  the  amount  of  the  kth  Zemike  polynomial 
Zk(m,  n )  present  in  the  optics.  These  coefficients  are  unique  to  the  imaging  system  being 
used.  To  create  an  accurate  model  of  the  PSF,  the  aberrations  present  in  the  telescope 
need  to  be  experimentally  measured.  Finding  the  values  for  these  coefficients  can  be 
difficult,  especially  in  a  three-mirror  optical  system,  but  Woods  has  described  a  method 
for  obtaining  the  Zemike  coefficients  values  [63].  Using  this  method,  the  SST  program  can 
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experimentally  estimate  the  first  11  coefficients.  These  values  can  change  over  time,  but 
average  values  are  used  in  this  research  to  simulate  the  aberrations  in  the  system.  Only 
the  first  1 1  coefficients  are  used  in  this  modei  for  two  reasons.  The  first  reason  is  that  a 
iarge  portion  of  the  power  in  the  aberrations  is  contained  in  the  iowest  order  poiynomiais. 
Secondly,  the  SST  system  is  undersampled,  and  any  high  frequency  distortions  will  not  be 
observable  in  th  sampled  images.  Table  3.1  shows  the  measured  coefficients  for  the  first  1 1 
poiynomiais. 


Table  3.1:  Average  measured  optical  aberration  Zernike  coefficients,  cij  for  the  SST. 


Coefficient 

Value  (Waves) 

Coefficient 

Value  (Waves) 

al 

2.09 

a7 

0.28 

Cl! 

-5.95 

a8 

-0.73 

a3 

-5.30 

a9 

0.36 

04 

6.89 

a  lo 

-0.48 

a5 

1.26 

an 

-0.16 

a6 

-0.28 

The  next  important  effect  to  modei  is  the  atmosphere.  The  atmosphere  is  well  modeled 
as  a  random  process  to  any  light  that  travels  through  it.  This  results  in  a  propagation  path 
for  the  light  that  has  varying  indices  of  refraction  in  space  and  time.  The  difference  in  path 
lengths  results  in  phase  front  distortions  and  the  image  is  not  formed  correctly  at  the  image 
plane.  A  telescope  collects  photons  over  a  specified  period  of  time  known  as  the  integration 
time,  T.  The  random  atmosphere  acts  differently  as  a  function  of  the  integration  time  used, 
ft  has  been  shown  that  over  long  integration  times,  T  »  100ms,  the  telescope  effectively 
averages  the  random  atmosphere  [17].  An  expression  for  this  average  random  atmosphere 
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or  long  exposure  Optical  Transfer  Function  (OTF),  HL,  is  derived  in  [17].  The  OTF  is  the 
Fourier  transform  of  the  PSF,  and  represents  the  frequency  response  of  the  system  to  an 
impulse. 


I 


HdfxJy )  =  exp 


-3.44 


ro 


(3.4) 


v  J  ) 

f x  and  fy  are  spatial  frequency  variables  at  the  image  plane,  A  is  the  center  wavelength 
of  the  telescope,  z  is  the  focal  length,  and  r0  is  the  seeing  parameter.  The  seeing  parameter 
is  a  variable  that  characterizes  the  intensity  of  fluctuations  in  the  atmosphere.  A  smaller 
r0  indicates  stronger  fluctuations,  implying  the  atmosphere  will  have  a  larger  effect  on  the 
telescope.  Typical  observed  ranges  for  r0  at  the  SST  site  in  the  White  Sands  Missile  Range 
are  5- 15cm.  Next,  the  Fourier  Transform  of  the  optical  PSF  and  long  exposure  OTF  are 
multiplied  in  the  frequency  domain  to  obtain  the  combined  OTF.  Then  the  inverse  Fourier 
is  performed  to  obtain  the  new  PSF  model,  h(x,y ): 


h(x,y)  =  T-1  (T[hopt(x,y)\HL(fx,fyj)  (3.5) 

The  model  described  up  to  this  point  is  accurate,  if  the  image  presented  to  the  detector 
is  sampled  completely.  In  order  to  sample  completely,  the  system  must  sample  the  image 
greater  than  the  Nyquist  frequency,  fs.  The  Nyquist  frequency  is  twice  the  cutoff  frequency, 
fc,  of  the  system,  and  is  the  rate  at  which  a  system  needs  to  be  sampled  to  prevent  any 
aliasing.  The  maximum  frequency  content  possible,  or  the  cutoff  frequency,  in  an  optical 
system  is  [34]: 


(3.6) 
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where  z  is  the  focal  length  of  the  system  and  d  is  the  diameter  of  the  aperture.  Using  the 
inverse  relationship  between  sampling  frequency  and  sample  spacing,  the  necessary  sample 
spacing  A.s  is  found  by: 


l  Az 

Av  =  —  =>  Av  <  — 
2 fc  ~  2d 


(3.7) 


If  an  optical  system  is  not  properly  sampled,  aliasing  occurs  and  data  is  lost.  A  value 
of  500nm  is  used  for  A  in  this  model  because  it  is  close  to  the  center  of  the  visible  range 
and  produces  an  integer  ratio  between  actual  pixel  size  and  Nyquist  sampling.  Solving 
equation  (3.7)  with  the  system  parameters  of  the  SST,  the  required  pixel  spacing  needed  to 
accurately  sample  the  intensity  at  the  CCD  array  is  Av  =  0.25/nn.  For  most  data  collections 
by  the  SST,  the  15//m  pixels  are  grouped  into  2x2  data  bins,  giving  an  effective  pixel  size 
of  A,  =  30/rm.  Since  Ar  >  As,  the  SST  is  not  sampled  at  Nyquist.  This  implies  that  a 
shift  in  where  the  image  is  formed  at  the  CCD  can  potentially  cause  differently  shaped 
PSFs.  To  be  able  to  simulate  this  effect,  shifting  is  performed  at  this  point  in  the  model, 
prior  to  adjustment  for  the  aliasing  caused  by  sampling.  The  effect  of  shifting  before  or 
after  modeling  the  sampling  effect  is  well-documented  in  [24].  The  author  demonstrates 
that  shifting  after  sampling  does  not  adequately  represent  the  effects  of  sub-pixel  location 
object  changes.  To  shift  the  PSF,  h(x,y),  two  variables  are  introduced;  a  distance  shift  in 
the  x  direction,  a,  and  y  direction,  co.  Since  the  PSF  is  created  in  a  large  zero  padded  matrix, 
utilizing  a  circular  shift  introduces  no  edge  effects. 


hs(x,y)  =  h(x  +  a,y  +  (o)  (3.8) 

After  the  shift  is  implemented,  the  difference  between  required  and  actual  sampling  is 
modeled  with  a  blurring  function.  Assuming  square  CCD  sensors,  this  blurring  function 
is  expressed  as  rectangles  in  the  x  and  y  dimensions  with  a  width  of  J3.  the  ratio  of  system 
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sampling  vs.  required  Nyquist  sampling.  This  more  accurate  model  for  the  PSF,  hd(x,y),  is 
produced  by  convolving  the  optical  PSF  from  equation  (3.8)  with  the  following  equation: 

hd(x,y)  =  ^  ^  hs(xx , yi)rect  j rect  yi  j  (3.9) 

xi  and  yi  are  convolution  variables  and  Nx  and  Ny  are  the  total  number  of  pixels  in  the 
V  and  y  direction  respectively,  while  hd(x,  y)  represents  the  convolved  and  downsampling 
effects  on  the  PSF  of  sampling  at  a  frequency  greater  than  Nyquist. 

This  PSF  model  is  used  for  this  simulation  because  it  combines  the  lens  aberrations 
with  the  actual  sampling  size  for  the  system,  further  increasing  the  fidelity  of  the  SST 
model.  Finally,  the  new  PSF  is  normalized  to  sum  to  one.  This  allows  the  PSF  to  be  scaled 
to  match  a  specified  object  intensity. 


hn(x,  y) 


hd(x,y ) 

NX  Ny 

X  X  hd(x,y ) 

X=1  X=1 


(3.10) 


hn(x,y)  is  the  final  modeled  PSF  used  in  this  chapter.  Depending  on  the  a  and  a>  shifts 
used,  the  values  and  distribution  of  hn(x,y)  will  change.  Next,  the  theory  and  motivations 
for  the  detection  algorithm  being  developed  is  described. 


3.2  Theory  Development 

In  this  section,  an  algorithm  is  developed  for  a  Multiple  Hypothesis  Test  (MHT). 
There  are  several  benefits  to  using  a  MHT,  along  with  a  strong  physical  motivation.  Using 
a  MHT  algorithm  presents  unique  computations,  which  are  addressed  in  this  section. 

3.2.1  Motivation  and  Selection  of  Alternate  Hypotheses. 

In  a  MHT,  M  different  potential  hypotheses  are  considered.  Each  hypothesis 
corresponds  to  a  particular  set  of  input  conditions.  In  this  case,  each  input  condition 
represents  a  modeled  PSF  viewed  through  the  SST  system.  As  described  in  section  3.1, 
sub-pixel  shifts  of  an  object’s  location  on  the  CCD  array  will  change  not  just  the  position, 
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but  the  shape  of  the  resulting  image.  This  effect  due  to  aliasing  will  cause  the  PSF  to  have 
different  shapes  depending  on  the  location  of  the  object  at  the  CCD  [24].  The  goal  of 
using  a  MHT  is  to  utilize  several  different  potential  PSF  distributions  to  achieve  a  higher 
probability  of  detection  for  an  object  that  may  be  located  in  one  of  the  sub-pixel  locations. 
Table  3.2  shows  a  list  of  potential  hypotheses  used  in  this  model. 


Table  3.2:  List  of  hypotheses  considered  and  a  description  of  object  location  at  the  CCD. 


Hypothesis 

Description 

Ho 

No  Object  Present 

H , 

a i,  a>i  Shifts 

h2 

a2,  co2  Shifts 

Hm-\ 

ctm-i,  Shifts 

By  picking  a  and  to  to  be  specific  sets  of  spatial  shifts  in  x  and  y,  a  model  for  objects 
located  in  the  corners,  center,  and  sides  of  a  pixel  can  be  developed.  Choosing  the  number 
of  hypotheses  and  their  position  location  is  an  important  decision.  As  described  in  Chapter 
2,  only  two  hypotheses  are  currently  used  by  the  SST  program.  These  two  hypotheses 
are  that  the  object  is  present  or  not  present.  Through  a  Likelihood  Ratio  Test  (LRT), 
the  test  is  reduced  to  examining  individual  pixels  and  determining  their  SNR.  If  a  space 
object  is  located  in  the  corner  of  a  pixel,  photons  may  spread  into  adjacent  pixels,  and 
the  peak  intensity  of  the  pixel  of  interest  is  lowered.  This  effect  lowers  the  SNR  value  of 
the  pixel  being  tested,  decreasing  the  probability  of  detecting  that  object.  To  mitigate  this 
effect,  correlation  algorithms  are  proposed  to  compare  not  just  a  single  pixel,  but  the  entire 
expected  PSF  [35].  This  concept  is  expanded  to  include  multiple  alternative  hypotheses 
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consisting  of  shifted  aliased  PSFs  that  are  used  as  potential  alternative  hypotheses  [24]. 
Figure  3.1  shows  a  graphical  representation  of  an  object  formed  on  a  single  pixel  with 
three  different  M  values,  along  with  the  true  object  location.  The  small  dot  indicates  the 
true  location  of  the  object,  and  the  lines  indicate  the  separation  of  the  pixel  into  its  decision 
space  areas. 


(A)  (B)  (C) 


Figure  3. 1:  A  graphical  example  of  the  division  of  a  pixel  into  sub-pixel  location  for  varying 
M  where  the  star  in  each  figure  represents  the  modeled  location  of  the  object  and  the  dot  is 
the  true  location.  (A)  Simple  Binary  test,  M  -  2  (B)  Multiple  hypothesis  test  with  M  =  6, 
with  sub-pixel  hypotheses  at  the  four  corners  and  center  of  the  pixel  (C)  A  highly  sampled 
pixel  space  with  M  -  121  hypotheses. 


A  BHT  is  illustrated  in  (A).  In  this  case,  the  object  is  assumed  to  be  in  the  center 
of  the  pixel.  Given  that  the  object  is  actually  located  in  the  bottom  left  comer,  the  point 
detector  or  correlation  BHT  might  not  adequately  capture  the  resulting  PSF.  This  can 
manifest  as  a  lower  probability  of  detection.  In  (B),  a  MHT  where  M  -  6  hypotheses  is 
illustrated.  In  this  case,  a  PSF  is  created  that  would  model  an  object  in  the  bottom  corner 
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adequately.  This  algorithm  would  likely  have  a  greater  probability  of  detection  due  to 
the  strong  correlation  between  the  modeled  and  observed  PSF.  Finally  in  (C),  a  densely 
sampled  pixel  is  illustrated.  Similarly  to  (B),  this  choice  of  M  would  produce  a  PSF  that 
would  likely  correlate  well  with  the  observed  PSF.  The  small  change  in  sub-pixel  locations 
used  in  (B)  and  (C)  would  not  present  a  noticeable  change  in  the  PSF. 

At  a  certain  point,  the  addition  of  an  increasing  number  of  hypotheses  introduces 
error  into  the  decision  criteria.  Since  the  PSFs  are  created  through  simulation,  many 
resulting  shapes  and  intensity  distributions  can  be  created.  For  a  large  number  of  M, 
the  alternate  hypotheses  may  have  a  higher  chance  of  resembling  random  noise  present 
in  the  system.  This  effect  can  potentially  increase  the  false  alarm  rate  of  the  detector, 
demonstrating  the  trade-off  between  the  number  of  hypotheses  and  the  detection  capability. 
In  addition,  the  creation  and  hypothesis  testing  of  the  additional  hypotheses  presents 
additional  computation  cost. 

Theoretically,  M  can  be  increased  to  infinity,  creating  a  continuously  sampled  space. 
This  result  would  lead  to  an  estimation  instead  of  a  detection  approach.  Following  this 
logic,  an  alternative  approach  to  the  MHT  presented  in  this  section  could  be  implemented. 
In  this  method,  an  estimate  is  found  for  the  sub-pixel  location  of  a  potential  object  at  every 
pixel.  Then  a  model  of  the  PSF  due  to  an  object  at  that  location  is  created,  and  a  BHT 
is  performed  based  on  that  estimate.  This  BHT  could  resemble  the  currently  used  SST 
algorithm.  This  method  would  result  in  performing  an  additional  estimation  algorithm  at 
every  pixel  location  in  addition  to  the  detection  algorithm,  which  is  a  significant  drawback. 

Another  drawback  with  an  estimation  of  the  location  approach  is  that  the  user  might 
not  be  interested  in  the  sub-pixel  location,  and  only  want  to  use  that  information  to  increase 
the  detection  performance.  In  an  estimation  problem,  the  sub-pixel  location  would  be 
considered  a  nuisance  parameter.  There  are  methods  of  performing  estimation  without 
estimating  nuisance  parameters,  described  by  Kay  [64],  but  they  are  utilized  when  the  end 
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goal  of  the  problem  is  to  estimate  information  about  a  signal.  The  space  object  detection 
algorithm  described  here  could  utilize  estimation  for  nuisance  parameters,  but  is  still  a 
detection  problem  at  the  core,  meaning  the  end  goal  of  the  user  is  to  determine  if  an  object 
is  present  or  not. 

An  additional  reason  that  this  problem  is  framed  as  a  MHT  is  that  it  allows  for  a  very 
specific  set  of  choices  for  sub-pixel  locations.  This  set  is  selected  as  the  most  likely  to 
occur,  and  will  have  the  largest  impact  on  the  intensity  and  shape  of  the  resulting  PSF. 
The  MHT  can  be  tailored  to  a  specific  optical  system  and  mission  requirements.  In  the  next 
section,  the  theory  involved  with  creating  and  implementing  a  MHT  algorithm  is  described. 

3.2.2  Detection  Algorithm  Based  on  Bayes  Risk. 

In  this  algorithm,  M  different  sub-pixel  locations  are  considered.  In  both  [25,  65]  a 
Bayes  criterion  is  shown  for  a  multiple  hypothesis  test,  including  methods  for  reducing  the 
test  to  a  usable  form.  A  Bayes  risk  hypothesis  test  details  how  to  select  the  hypothesis  that 
results  in  the  minimal  risk,  R,  on  average.  The  hypothesis  that  results  in  the  smallest  risk 
is  the  most  likely  to  match  the  true  hypothesis. 

M- 1  M-l  r 

nkCik  I  p(D\Hk)dD  (3.11) 

i=0  k= o  J Zi 

In  the  risk  expression,  nk  is  the  a  priori  probability  for  the  kth  hypothesis,  or  how 
likely  it  is  to  occur.  Cik  is  the  cost  associated  with  choosing  the  ith  hypothesis  when  the  k,h 
hypothesis  has  occurred.  The  cost  is  a  number  between  zero  and  one  that  indicates  how 
big  of  an  error  is  made,  with  zero  being  no  error  and  one  being  the  most  significant  type 
of  error.  The  conditional  PDF  of  the  received  data  given  Hk  occurred  is  p(D\Hk).  Zk  is 
observation  space  for  the  k,h  hypothesis. 

To  reduce  this  expression  into  a  form  that  can  be  implemented,  the  null  hypothesis 
space  Z0  is  separated  into  its  non-overlapping  regions  Z0  -  Z  -  Zx  -  ...  -  ZM_  | .  Next, 
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noting  that  integrals  over  the  entire  decision  space  Z  equal  one,  and  grouping  similar  terms 
of  Zk  gives  the  following  relation: 


m— 1  m— 1  rM- 1 

k = Yj  ^ckk +Yj  I  Tjnk  (Cik  - Ckk)  p(D lHk)dD-  (3-12) 

k  i  v  fc=o 

z'  k+i 

The  first  summation  term  in  equation  (3.12)  is  a  constant  risk,  implying  that  it  is 
independent  of  the  choice  of  region  i,  and  can  be  excluded  from  the  decision  criteria.  The 
cost  for  selecting  a  correct  hypothesis,  Ckk,  is  assumed  to  have  no  associated  risk  and  is 
defined  to  be  zero.  The  minimum  Bayes  risk  is  then  found  by  choosing  the  region  i  that 
results  in  the  smallest  risk  K.  The  i  that  results  in  the  minimum  risk  is  the  hypothesis  H, 
which  minimizes  Bayes  risk.  Determining  this  region  i  is  implemented  by  selecting  the 
integrand  over  Z,  with  the  smallest  value.  This  results  in  the  following  decision  criteria: 


(m-  i  'I 


Hi  =  min 

ieO:M-l 


nkCikp(D\Hk) 


(3.13) 


k= 0 

V  k+i  ) 

By  using  equation  (3.13)  directly,  a  detection  algorithm  can  be  implemented.  If  H0 
produces  the  minimum  risk  value,  then  the  algorithm  says  no  object  is  present.  On  the 
other  hand,  if  any  other  hypothesis  produces  the  minimum  risk,  it  is  decided  that  an  object 
is  present  at  the  sub-pixel  location  hypothesis  providing  the  minimum  risk  value.  This 
MHT  can  be  reduced  further  through  specific  choices  for  costs  and  priors,  or  implemented 
directly.  If  it  is  implemented  directly,  the  MHT  is  flexible,  because  any  combination  of 
distributions  for  p(D\Hk),  costs,  or  priors  can  be  used  without  changing  the  basic  structure 
of  the  algorithm.  A  drawback  of  this  method  is  that  finding  the  smallest  hypothesis  requires 
non-linear  functions  and  is  computationally  complex.  In  this  research  effort,  both  a  direct 
implementation  and  a  reduced  algorithm  based  on  specific  choices  are  investigated.  To 
further  reduce  the  MHT  from  equation  (3.13),  a  few  additional  details  are  needed,  including 
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a  model  for  the  conditional  probabilities  p(T)\Hk),  and  a  structure  for  the  prior  probabilities 
nk  and  costs  Cik.  These  factors  are  investigated  in  the  next  sections. 

3.2.3  Model  for  Received  Data  Given  Hypothesis  Hk. 

An  important  component  for  a  hypothesis  test  is  having  a  model  for  the  received  data. 
It  is  well-known  that  photon  arrival  times  of  light  is  a  Poisson  process  [17].  Photon  arrival 
times  are  not  the  only  source  of  noise  present  in  the  SST  system.  The  CCD  can  introduce 
several  different  sources  of  error,  including  dark  current  and  spatial  non- uniformity  across 
pixels  [58,  66].  Under  the  H0  hypothesis,  when  no  object  is  present,  the  measured  data 
contains  only  low  intensity  background  light  and  noise  from  the  CCD.  Since  the  received 
data  is  dominated  by  CCD  noise,  it  is  assumed  to  be  Gaussian.  In  the  other  M-l  hypotheses, 
there  is  considered  to  be  an  object  present,  and  therefore  Poisson  noise  will  contribute  more 
significantly  towards  the  overall  noise  level.  Despite  the  additional  Poisson  noise,  it  is 
assumed  that  the  Gaussian  noise  sources  still  dominate  in  these  alternate  hypotheses.  As 
a  result,  the  received  data  given  these  hypotheses  is  also  considered  to  be  Gaussian.  As 
described  in  Chapter  2,  several  research  efforts  have  investigated  potential  alternate  PDFs 
for  these  hypotheses  [22,  35].  Making  these  stated  assumptions,  the  conditional  PDF  for  a 
single  pixel  of  received  data  given  Hk  is  expressed  as: 


1  /  (D  -  I)2  \ 

P{Dm=  vs7exp(-— )  <3'14) 

/  is  the  expected  mean  of  the  received  data.  If  no  object  is  present,  H0,  the  expected 
mean  is  only  the  background  photons  observed,  I  =  B.  In  the  other  hypotheses,  /  is  the 
object  intensity  6  plus  the  constant  background,  I  =  6  +  B.  In  both  cases,  cr  is  the  standard 
deviation  of  the  received  data.  Assuming  each  pixel  in  the  image  is  independent,  the  joint 
conditional  probability  of  received  image  data  D  given  the  kth  hypothesis  is: 


NyNy 


P(D\Hk)  -  (2ncr2^  2  exp 


t  (P(x ,  y)  -  I(x,  y))2 

V  X=l  y=l 


2cr2 


(3.15) 
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In  this  case,  I(x, y)  =  B  for  H0  and  I(x, y)  =  Ohkn{x,y)  +  B  for  Hk.  Additionally,  hkn(x,y ) 
is  the  normalized  model  for  the  PSF  of  the  kth  hypothesis. 

3.2.4  Selecting  a  Cost  Structure. 

To  implement  the  MHT,  it  is  necessary  to  choose  values  for  the  costs,  Cjk,  and 
priors,  nk.  Often  if  there  is  not  strong  evidence  otherwise,  an  equal  cost  and  uniform 
prior  probability  scheme  is  chosen.  This  cost  scheme,  a  Equal-Cost  Equal-Prior  (ECEP) 
assignment,  weighs  all  errors  with  the  same  amount  of  cost.  For  example,  a  false  alarm  has 
the  same  cost  as  an  error  in  sub-pixel  detection  and  a  missed  detection.  It  also  assumes  that 
all  hypotheses  occur  with  the  same  probability. 


I  0,  i  =  k  1 

Co,  =  \  ,  nk  =  —  (3.16) 

[  1,  i  *  k  M 

In  certain  situations,  this  cost  structure  accurately  represents  the  true  state  of  the 
system.  In  a  SDA  environment,  it  is  much  more  probable  that  there  is  no  object  present, 
Hq,  than  having  an  object  present,  Hke i:M_i  due  to  the  relatively  small  number  of  detectable 
objects  in  a  large  search  volume.  As  described  in  Chapter  2,  a  detection  algorithm  utilizing 
uniform  priors  and  equal  costs  is  implemented  [24].  One  reason  for  its  use  is  its  ability  to 
reduce  the  algorithm  to  a  sufficient  statistic  or  SNR  test.  Using  SNR  is  a  common  practice 
in  other  SDA  efforts  [21].  This  approach  gives  the  author  the  ability  to  directly  compare 
detection  algorithm  performance  to  the  point  detector  used  by  the  SST.  Using  these  equal 
cost  and  prior  assumptions  equation  (3.13)  reduces  to  the  following  LRT  comparison: 


Le  =  max 
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(3.17) 
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B  is  the  estimated  background  intensity  and  is  found  by  determining  the  median  value 
of  all  the  pixels  in  the  15x15  window,  cr  is  the  standard  deviation  of  the  received  data.  To 
calculate  cr,  the  standard  deviation  of  the  15x15  window  is  found.  Next,  the  outliers  with 
intensity  more  than  3  standard  deviations  above  the  estimated  background  are  removed. 
A  new  cr  is  calculated  with  these  outliers  removed.  More  details  about  this  process  are 
described  in  Section  3.4.  crk  is  a  normalizing  factor  that  is  included  to  ensure  that  X,, 
has  the  same  standard  deviation  as  the  original  point  detector.  The  value  for  crk  is  found 
by  summing  the  squared  PSF  of  the  klh  hypothesis.  The  SNR  for  each  hypothesis  k  is 
computed,  and  the  maximum  SNR  is  compared  to  a  threshold  rM_ary. 

Alternatively,  this  model  implements  a  cost  scheme  with  Unequal-Cost  Equal-Prior 
(UCEP)  assumptions,  where  errors  of  different  types  are  weighted  differently.  A  third  type 
of  error  besides  false  alarm  and  missed  detection  is  introduced.  This  third  error  occurs 
when  an  object  is  detected,  but  is  misclassified  as  having  the  incorrect  sub-pixel  location. 
There  is  no  cost  associated  with  this  “mixed-detection”  event.  This  cost  is  useful  in  a 
situation  where  the  detection  of  the  object  is  more  highly  valued  than  a  precise  location,  as 
is  the  situation  for  the  SST.  For  example,  in  the  SST  program,  the  detection  results  feed 
into  other  data  processing  steps  that  perform  clustering  and  locating.  Missed  detections  and 
false  alarms  are  still  considered  to  have  a  cost  of  one.  To  compare  to  previous  algorithms, 
an  equal  prior  probability  for  nk  assumption  is  made  for  this  MHT. 
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(3.18) 


Implementing  the  costs  from  equation  (3.18)  into  the  Bayes  risk  from  equation  (3.13) 
and  reducing  gives  the  following  detection  criteria: 
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vfc'pCDI^o): 

Using  the  assumed  Gaussian  distributions  described  in  equation  (3.15)  for  /?(D| ///.),  the 
ratio  of  conditional  densities  can  be  reduced.  Combining  the  exponentials  and  canceling 
like  terms  gives  the  following  relation  for  J2U: 
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(3.20) 


Unlike  an  equal  cost  approach,  this  MHT  does  not  simplify  to  a  SNR  calculation 
due  to  the  summation  of  all  of  the  potential  alternate  hypotheses  which  are  exponentials. 
This  demonstrates  both  how  additional  information  is  being  used  to  make  the  detection 
decision,  and  how  the  cost  structure  reduces  the  test  to  a  decision  between  only  Hi  and  H{). 
This  implies  that  no  sub-pixel  information  is  inferred  from  this  test.  Another  distinguishing 
factor  between  an  equal  and  unequal  cost  approach  is  that  the  intensity  of  the  object  9k  is 
now  needed  to  determine  JLU.  The  intensity  6k  is  the  object  intensity  used  in  the  current 
MHT  method  presented  in  [24].  There  are  different  intensities  for  each  hypothesis  because 
the  shifting  results  in  more  or  less  spreading  of  the  PSF,  causing  intensity  changes.  In  the 
previous  equal  cost  MHT  algorithm,  6  is  not  defined  directly,  but  is  chosen  by  the  selection 
of  a  threshold  TM-ary  since  9  is  included  in  the  threshold. 
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(3.21) 


The  algorithm  presented  in  equation  (3.20)  is  a  reduced  form  of  the  decision  criteria 
shown  in  equation  (3.13),  but  is  specific  to  a  set  of  costs.  To  evaluate  the  performance  of 
equal  and  unequal  costs  in  the  detection  of  objects,  the  non-reduced  equal  cost  approach 
presented  in  (3.13)  is  compared  to  the  reduced  unequal  cost  algorithm  derived  in  (3.20). 
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There  are  benefits  and  drawbacks  to  using  both  the  ECEP  and  ETCEP  algorithms. 
ETsing  the  previously  developed  equal  cost  algorithm  shown  in  equation  (3.17)  gives  both 
detection  capability  and  sub-pixel  location  information.  The  largest  hypothesis  SNR  is 
compared  to  the  threshold,  and  if  an  object  exceeds  that  threshold  it  is  considered  to  be 
in  that  sub-pixel  location.  The  additional  position  information  can  be  useful  in  improving 
object  tracking.  A  drawback  of  this  algorithm  is  that  each  SNR  calculation  is  only  based 
on  one  PSF  model.  Each  of  the  M  hypotheses  are  correlated  against  the  received  data,  D, 
and  are  used  to  choose  the  maximum.  This  equal  cost  method  of  detection  is  a  series  of 
binary  hypothesis  tests  fused  together.  The  alternative  ECEP  algorithm  presented  with  this 
research  is  a  direct  implementation  of  equation  (3.13),  resulting  in  similar  benefits  to  the 
previously  developed  approach.  The  benefit  to  the  new  implementation  is  that  it  is  flexible 
in  its  ability  to  change  the  prior  probabilities  or  costs  if  desired. 

In  contrast  to  the  ECEP,  the  UCEP  algorithm  consists  of  a  combination  of  all  potential 
hypotheses  for  the  selection  of  k.  This  implies  the  data  from  each  hypothesis  is  combined  to 
make  a  single  decision.  A  drawback  of  the  unequal  cost  detection  algorithm  from  equation 
(3.20)  is  that  it  is  unable  to  be  reduced  to  a  sufficient  statistic  such  as  SNR.  As  a  result,  there 
are  additional  computation  complexities  which  may  limit  performance  or  implementation 
operationally. 

A  significant  difference  with  this  new  approach  is  the  space  in  which  decisions  are 
made.  Previous  work  has  reduced  the  hypothesis  test  to  a  sufficient  statistic,  namely  SNR. 
Using  the  SNR,  it  is  more  straightforward  to  determine  the  probability  of  detection.  On  the 
other  hand,  using  an  UCEP  algorithm  (3.20)  gives  a  complicated  relationship  between  M 
hypotheses  in  risk  space.  As  part  of  this  effort,  a  method  to  determine  Pd,  P /,  and  create 
ROC  curves  in  risk  space  is  developed.  These  derivations  are  included  in  section  3.4.  Next, 
the  experiment  used  to  collect  data  is  described. 
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3.3  Experiment  Description 

To  measure  the  performance  of  the  detection  algorithms,  data  collected  from  the 
SST  system  is  analyzed.  In  operational  use,  the  telescope  will  observe  varying  intensity 
objects.  To  demonstrate  that  a  new  detection  algorithm  performs  better  than  the  currently 
implemented  BHT,  objects  of  high,  intermediate,  and  low  intensity  need  to  be  investigated. 
An  improvement  in  detection  of  any  of  these  types  of  objects  will  enhance  the  overall 
capabilities  of  the  SST  system. 

One  approach  to  test  the  MHT  would  be  to  find  a  variety  of  known  and  cataloged 
space  objects  of  varying  intensity,  and  image  these  objects  with  the  SST.  It  may  be  difficult 
to  find  enough  objects  of  particular  intensity  without  a  large  search.  It  could  also  require 
a  long  period  of  search,  or  multiple  nights  to  capture  these  objects.  During  that  time, 
important  parameters  may  change,  including  the  seeing  r0  or  optical  aberrations  in  the 
telescope.  Instead,  this  research  uses  a  data  set  of  the  SST  observing  the  Geosynchronous 
Earth  Orbit  (GEO)  communications  satellite  ANIK-F1  as  it  enters  into  an  eclipse.  The  SST 
is  programmed  to  track  the  orbital  elements  of  the  satellite  as  it  enters  an  eclipse  caused 
by  the  earth.  Collecting  images  in  this  manner  causes  the  satellite  to  change  intensity 
as  a  function  of  time,  transitioning  from  bright  to  dim.  Collecting  images  of  this  event 
reduces  the  requirement  of  searching  several  areas  of  sky  and  locating  objects  of  specific 
magnitudes. 

The  data  set  described  here  was  originally  collected  to  observe  the  performance  of 
another  detection  algorithm  presented  in  [24].  The  result  of  this  observation  is  a  set  of  data 
including  a  space  object  that  transitions  from  easily  detectable  to  very  difficult  to  detect 
over  the  course  of  a  few  minutes,  or  a  few  hundred  images.  This  experiment  was  conducted 
during  the  vernal  equinox  in  2012  and  was  repeated  over  several  nights  from  28  Feb  2012 
to  23  Mar  2012. 
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3.3.1  Locating  and  Selecting  False  Alarm  Data. 

To  completely  compare  a  set  of  algorithms,  the  detection  performance  should  be  noted 
for  different  false  alarm  rates.  Pf  is  the  probability  that  when  the  algorithm  is  presented 
with  data  containing  no  object,  it  indicates  that  an  object  is  present.  To  find  a  value  for  Pf, 
captured  data  containing  no  object  is  needed.  There  are  two  potential  methods  to  create 
data  with  no  object  present. 

One  method  is  to  generate  independent  Poisson  random  variables  at  each  pixel  with 
a  mean  of  the  background  intensity.  This  simulated  data  would  give  an  approximation  of 
actual  data  that  might  be  observed  by  the  SST  when  no  object  is  present.  Alternatively,  this 
research  uses  collected  SST  data  to  be  as  realistic  as  possible.  Due  to  the  potential  for  dim 
objects  at  any  point,  it  can  be  difficult  to  determine  if  there  is  truly  an  object  such  as  a  star 
in  any  subset  of  data.  There  may  be  an  object  present,  but  barely  above  the  noise  floor  and 
not  noticeable  to  the  human  eye.  To  ensure  a  set  of  data  is  used  with  no  object,  a  patch 
of  sky  is  collected  and  averaged  over  multiple  frames.  Since  the  SST  system  is  tracking  a 
satellite  in  GEO,  most  of  the  objects  not  in  GEO  will  move  at  the  sidereal  rate.  This  rate 
translates  to  8  pixels  of  motion  per  frame  or  greater. 

To  generate  a  set  of  false  alarm  data,  a  15x15  pixel  patch  of  observed  sky  is  tracked 
through  60  consecutive  frames.  The  60  frames  are  then  averaged  to  give  a  longer  effective 
integration  time.  If  there  were  a  dim  celestial  object  consistently  present  in  this  portion  of 
sky,  it  would  become  more  intense  as  the  noise  is  reduced  through  the  averaging  process. 
Through  several  trials,  a  portion  of  sky  is  found  and  documented  with  no  pixel  outside  of  3 
standard  deviations  of  the  mean  in  this  averaged  data  set.  This  implies  that  there  is  a  high 
likelihood  that  all  of  the  fluctuations  are  due  to  noise,  and  no  object  is  present  in  the  data 
set.  This  data  set  is  then  used  to  calculate  an  upper  bound  of  the  false  alarm  rates  used  in 
the  results. 
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3.4  Results  and  Analysis 

Each  image  collected  by  SST  consists  of  6144x4096  pixels  that  are  2x2  binned.  To 
perform  the  detection  algorithm  in  a  tractable  amount  of  time,  a  small  subset  of  the  image 
around  the  satellite  of  15x15  pixels  is  extracted  from  each  frame.  This  allows  for  the 
algorithm  to  quickly  process  a  MHT  on  the  center  pixel,  but  still  includes  enough  pixels  to 
capture  the  entire  PSF  and  a  sufficient  amount  of  background  pixels  to  calculate  accurate 
background  statistics.  In  operational  use,  this  15x15  window  surrounding  the  pixel  of 
interest  would  slide  through  the  image  as  each  new  pixel  is  tested.  For  this  analysis,  the 
window  surrounding  the  center  pixel  of  ANIK-Fl’s  PSF  is  investigated  in  each  frame. 

Another  important  factor  of  the  data  analysis  is  outlier  removal.  Within  each  15x15 
window  there  may  be  objects  in  addition  to  the  object  of  interest.  These  “nuisance  objects” 
are  detrimental  to  the  test  being  performed  on  the  center  pixel.  When  the  background,  B ,  is 
estimated  by  finding  the  median,  the  nuisance  objects  inflate  the  value  of  the  background 
statistics,  making  the  object  of  interest  seem  more  dim  relative  to  the  background.  This 
effect  causes  a  decrease  in  detection  performance. 

These  objects  are  removed  in  a  two-step  process  similar  to  the  method  described  in 
[24].  First,  the  background  and  variance  is  estimated  for  the  entire  15x15  window.  Next, 
the  individual  pixels  are  compared  against  the  estimated  background.  Any  pixels  that  are  3 
standard  deviations  or  greater  than  the  background  are  removed.  This  does  not  include  the 
PSF  of  the  object  of  interest.  Finally,  the  standard  deviation  is  calculated. 

3.4.1  Detection  and  False  Alarm  Computations. 

As  mentioned  in  section  3.2.4,  a  unique  method  of  computing  Pci  and  Pj  is  proposed 
in  this  section.  Fooking  at  equation  (3.13),  there  are  two  important  cases  in  determining 
detections  and  false  alarms.  These  cases  are  the  value  of  the  algorithm  when  H0  is  chosen, 
tj.  and  the  minimum  value  of  the  rest  of  the  M- 1  other  hypothesis,  (p ■  To  give  a  similar 
comparison  to  the  UCEP  algorithm,  the  logarithm  is  taken  of  both  cases  to  move  the 
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variables  into  a  log  space.  The  data  being  processed,  D,  is  contained  in  the  exponential 
function.  As  a  result  it  is  common  in  detection  and  estimation  to  make  decisions  in  log 
space.  The  expression  for  these  random  variables  are  defined  by  the  two  expressions  below: 
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The  subscript  p  signifies  which  case  the  variables  were  generated  under.  These  random 
variables  are  generated  by  two  different  conditions.  The  first  condition  p  =  1  is  when  an 
object  is  present  in  the  data.  In  this  case,  D  equals  the  15x15  frame  of  data  surrounding 
ANIK-F1.  This  data  is  used  to  generate  the  probability  of  detection.  The  sub-pixel  location 
can  also  be  determined  by  noting  which  i  provides  the  minimum  from  equation  (3.23).  The 
accuracy  or  benefits  of  using  this  estimate  for  the  sub-pixel  location  of  the  object  is  not 
investigated  in  this  research.  The  second  condition,  p  -  0,  is  when  image  data  containing 
no  objects,  as  described  in  section  3.3.1,  is  used  for  D.  Processing  this  data  gives  the  false 
alarm  rate  of  the  test.  To  create  a  single  frame  detection  from  this  MHT,  the  two  variables 
can  be  compared  with  the  following  relation: 


Tip  ^  <Pp  +  t  (3.24) 

Ho 

If  (pp  is  less  than  r/p,  an  object  is  considered  to  be  present;  if  the  opposite  is  true,  no 
object  is  present.  The  threshold  r  can  be  changed  to  allow  the  user  to  achieve  a  desired  Pj. 
Implementing  equation  (3.24)  would  yield  one  instance  of  a  detection  decision  for  a  single 
frame.  This  value  may  not  be  representative  of  how  often  that  outcome  occurs.  Instead, 
a  sequence  of  frames  are  used  to  generate  a  probability  of  detection,  Pd,  at  that  particular 
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intensity.  Finding  the  probability  of  detection  relies  on  knowledge  of  the  distribution  of  the 
variables  <pp  and  r/p. 

As  described  in  section  3.2.3,  the  conditional  probabilities  p  (D|F4)  are  assumed  to  be 
Gaussian.  Investigating  the  values  of  the  variables  over  the  collected  data  shows  that  the 
random  variables  appear  to  follow  a  Gaussian  distribution.  To  verify  this,  a  sliding  window 
of  15  frames  is  used  to  compute  the  statistics  for  rj  and  0.  An  interesting  trade-off  exists 
due  to  the  object  continuously  losing  intensity.  Ideally,  to  ensure  the  best  fit  for  statistics, 
a  large  window  would  be  used.  This  would  allow  high  confidence  in  both  the  Gaussian 
distribution  assumption  and  the  calculated  mean  and  standard  deviations  of  the  data.  A 
drawback  of  using  too  large  of  a  window  is  that  the  intensity  of  the  satellite  decreases 
too  much  between  the  beginning  and  the  end  of  window.  This  results  in  data  that  appears 
Gaussian,  but  has  a  non-constant  mean.  In  this  case,  the  mean  value  of  rjp  and  <pp  is  a 
function  of  the  frame  number,  or  object  intensity.  To  avoid  this,  a  sliding  window  of  15 
frames  is  used.  This  window  was  large  enough  to  generate  accurate  statistics,  but  not  long 
enough  to  cause  the  intensity  of  ANIK-F1  to  change  significantly.  Making  this  assumption, 
the  mean  and  variance  statistics  can  describe  the  behavior  of  both  rjp  and  (pp.  Figure  3.2 
shows  a  plot  of  Tip  and  <pp  values  when  p  -  1  over  a  15  frame  window. 

For  the  case  shown  in  Figure  3.2,  the  value  for  <f>i  is  always  less  than  the  value  of  rp. 
If  a  single  point  detection  decision  is  used,  a  detection  would  be  made  at  every  point, 
assuming  a  threshold  value  of  r  =  1.  Designating  a  value  of  1  for  a  detection  or  0 
for  no  detection  at  each  frame  could  give  one  instance  of  the  detection  capability  of  the 
algorithm.  Using  a  single  point  approach  as  described  may  not  represent  the  typical  or 
average  detection  performance  given  slightly  different  noise  or  other  fluctuations.  Instead, 
the  means  p(/,  and  pn,  and  standard  deviations  (r(j}  and  cr;;  are  used  to  determine  how  probable 
a  detection  is  in  general.  To  generalize  the  decision  rule  in  equation  (3.24)  to  more  than 
one  instance,  the  probability  of  detection  or  false  alarm  Pe  is  created.  Pe  represents  the 
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Figure  3.2:  Equal  cost  algorithm  variables  with  an  object  present,  77 1  and  (p\.  These 
variables  are  observed  over  a  15-frame  window,  and  the  means  and  standard  deviations 
are  used  to  determine  Pc\  at  frame  800.  This  data  was  collected  on  13  Mar  2012. 


probability  that  the  random  variable  (pp  is  less  than  the  random  variable  r]p.  In  general,  the 
probability  that  one  random  variable  is  greater  than  another  is  not  directly  computable,  so 
instead  an  intermediate  conditional  probability  for  a  given  r/p  is  calculated. 

vP 

Pe(rjp)  =  Pd,((pp  <  pp\pp)  =  J' p </,((/) p)d(f)p  =  F^pp)  (3.25) 

—00 

The  conditional  probability  p^cpp  <  PP\pP)  causes  the  overall  probability  rP(,  to  become 
a  function  of  pk-  The  probability  that  a  random  variable  is  any  value  below  a  threshold  is 
defined  as  a  Cumulative  Distribution  Function  (CDF).  The  CDF  F^pp)  is  calculated  by 
integrating  the  PDF  for  <p,  pUpp),  from  -00  to  the  threshold  pp.  Evaluating  Ft/)(pp)  gives  the 
probability  of  detection  or  false  alarm  for  a  single  given  value  of  pp.  To  determine  rP(,  not 
as  a  function  of  pp,  Bayes  Theorem  is  utilized,  and  the  distribution  is  integrated  over  all 
possible  values  of  pp. 
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p,,(r]p)  is  the  PDF  for  the  random  variable  77  evaluated  at  tjp.  At  this  point  an  analytical 
solution  for  Pe  is  not  easily  determined.  To  calculate  Pe.  the  integral  over  r\p  is  solved 
numerically  using  a  Riemann  sum. 


Pe  =  ^  F p)pv(7ip)Ar]p  (3.27) 

ip 

To  analytically  solve  the  integral,  a  finite  range  of  t]p  values  are  needed.  i]p  is  varied 
from  -200  to  200  with  A///;  =  1.  The  range  of  rjp  is  selected  to  be  large  enough  to  cover 
potential  function  values.  A  similar  approach  is  used  to  determine  Pd  and  Pf  for  the  unequal 
costs  algorithm  X„  in  equation  (3.20).  One  significant  difference  is  that  there  are  no  longer 
two  random  variables  to  compare.  Instead  there  is  one  random  variable  and  a  constant 
threshold.  As  a  result,  the  probabilities  Vu  can  be  calculated  by  comparing  to  a  fixed 
threshold  T„.  Figure  3.3  is  a  plot  of  the  value  of  Lu  calculated  as  a  function  of  frame 
number.  The  value  for  for  this  window  is  always  above  the  threshold,  assuming  T  =  0. 

The  statistics  of  a  sliding  15  frame  window  are  calculated  for  the  mean  /U£  and  standard 
deviation  <xr .  Pu  is  then  found.  The  probability  of  false  alarm  or  detection  is  defined  as  the 
probability  that  the  random  variable  is  greater  than  the  threshold  r„. 


r„ 

Pu  =  l-  F£u( r„)  (3.28) 

Utilizing  the  fact  that  the  PDF  integrates  to  one  over  its  entire  probability  space,  the 
probability  Vu  can  be  found  by  subtracting  the  CDF  of  X«  from  one.  Up  until  this  point, 
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Figure  3.3:  Unequal  algorithm  £u  for  a  15  frame  window  centered  at  800  for  13  Mar  2012. 
The  mean  and  standard  deviation  are  calculated  to  give  the  probability  of  detection  or  false 
alarm  Pu. 


the  probabilities  Ve  and  Vu  have  been  left  generalized.  The  benefit  is  that  the  calculations 
for  Pd  and  Pf  are  very  similar  and  only  the  input  conditions  are  varied.  Since  Pd  is  the 
probability  of  choosing  0  when  (/>  is  true,  and  Pf  is  the  probability  of  choosing  (/>  when 
T]  is  true,  the  expression  for  Ve  and  Vu  can  be  used  for  both  detection  and  false  alarm 
probabilities.  By  varying  p,  or  whether  an  object  is  present  or  not,  both  probabilities  can 
be  calculated. 


[  Pd ,  P  =  1 

Pe,u  =  (3-29) 

1^/.  F  =  ° 

Using  equations  (3.27)  and  (3.28),  values  of  Pd  and  Pf  based  on  the  algorithms 
described  in  section  3.2.4  are  calculated  from  the  collected  data.  The  first  result  investigated 
is  the  probability  of  detection  at  the  false  alarm  rate  currently  used  by  the  SST  system.  The 
first  step  is  to  verify  the  results  against  previous  methods.  After  this  has  been  verified, 
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conclusions  can  be  made  about  any  changes  or  improvements.  This  comparison  is  between 
an  equal  cost  MHT  approach  originally  presented  in  [24],  an  unequal  cost  MHT,  and 
the  currently  used  binary  point  detector.  Using  an  equal  cost  MHT  but  calculating  P d 
differently  and  achieving  similar  results  for  a  given  false  alarm  rate  with  the  same  collected 
data  would  verify  this  new  method. 


(A)  (B)  (C) 


Figure  3.4:  Pd  curves  for  ANIK-F1  as  it  enters  eclipse  on  3  consecutive  nights  (A)  13  Mar 
2012  (B)  14  Mar  2012  (C)  15  Mar  2012.  The  solid  black  line  is  the  UCEP  algorithm,  the 
blue  dashed  line  is  the  ECEP  algorithm,  the  red  dotted  and  dashed  line  is  a  binary  matched 
filter,  and  the  brown  dotted  line  is  a  point  detector.  These  detection  curves  are  generated 
with  Pf  =  4.56e  -  10  and  M  =  10. 


Figure  3.4  shows  the  probability  of  detection  curves  for  data  collected  from  the  SST 
on  three  consecutive  nights.  The  first  night  (A)  was  collected  13  Mar  2012.  The  UCEP 
algorithm  outperforms  the  ECEP  algorithm  across  all  of  the  frames,  with  a  difference 
range  of  10-60  percent  depending  on  the  frame  number.  The  UCEP  algorithm  also  shows 
increases  of  40-90  percent  over  the  currently  implemented  point  detector.  Both  the  ECEP 
and  UCEP  algorithms  start  close  to  Pd  =  1,  decrease  quickly  as  the  satellite  enters  the 
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eclipse,  but  the  UCEP  algorithm  only  drops  to  approximately  30  percent  detection  over 
the  same  period  of  the  eclipse.  From  [24],  the  ECEP  MHT  crosses  50  percent  detection  at 
approximately  frame  840,  while  the  new  equal  cost  MHT  developed  crosses  at  frame  825. 
The  point  detector  crosses  at  approximately  frame  780  and  the  UCEP  algorithm  does  not 
cross  50  percent  until  frame  875.  Similar  trends  are  present  in  the  other  nights,  with  UCEP 
MHT  Pd  improvements  between  5-65  percent  in  (B)  and  5-70  percent  in  (C)  over  the  ECEP 
MHT  and  10-90  percent  in  (B)  and  10-95  percent  in  (C)  over  the  point  detector. 

In  all  three  days  processed,  the  object  goes  from  easily  detectable  to  low  Pd  across 
approximately  100  frames.  In  [24]  the  transition  happens  across  approximately  85  frames. 
The  new  ECEP  MHT  seems  to  match  £e  from  equation  (3.17),  the  method  and  results 
presented  in  [24]  with  a  small  number  of  differences.  One  of  these  differences  is  that 
the  preprocessing  of  the  data  before  implementation  into  the  algorithm  may  be  different. 
Another  contributing  factor  is  whether  stars  in  close  proximity  to  ANIK-F1  are  removed 
or  not.  Stars  entering  the  small  frame  being  processed  can  negatively  impact  the  detection 
performance  by  driving  up  the  noise  statistics  for  the  window.  Another  difference  is  that 
the  data  is  processed  with  different  methods.  In  this  new  method,  a  SNR  statistic  is  not 
used  and  the  data  is  processed  in  risk  space.  As  a  result,  the  random  numbers  are  combined 
in  a  different  and  non-linear  method,  giving  a  slightly  different  performance. 

Another  difference  is  the  method  of  computing  Pf.  Actual  collected  data  selected 
to  have  no  object  present  is  used  instead  of  simulating  a  false  alarm  rate.  Using  the 
collected  false  alarm  data  described  in  section  3.3.1  with  the  algorithm  to  generate  false 
alarm  probability  as  opposed  to  estimating  Pf  from  the  threshold  gives  a  more  accurate 
measure  of  the  true  false  alarm  rate  of  the  tests.  Pf  is  calculated  to  be  7.3e-4  for  the  equal 
cost  algorithm  using  the  estimated  0,  based  on  the  threshold  of  Tm.^  =  6.22  from  [24]. 
The  false  alarm  rates  generated  by  running  the  ECEP  and  UCEP  MHT  are  not  equal.  The 
calculated  Pf  for  the  unequal  MHT  is  found  to  be  2.8e-8.  To  give  an  even  comparison 
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of  detection  across  frames,  a  threshold  r  value  of  -27.5  is  used  in  the  UCEP  algorithm, 
to  match  the  false  alarm  rates  between  the  two  algorithms.  Once  the  algorithms  have  the 
same  Pf,  the  three  consecutive  frame  rule  is  taken  into  effect.  As  mentioned  previously,  the 
SST  requires  the  current  BHT  point  detector  to  observe  a  detection  over  three  consecutive 
frames  to  consider  an  object  present.  This  helps  reduce  the  false  alarm  rate  due  to  single 
image  noise  spikes  or  other  non-constant  sources.  The  measured  Pj  and  Pd  rates  are  cubed 
to  account  for  this  effect. 

To  quantify  the  percent  improvement  in  detection  performance  with  an  additional 
metric,  the  apparent  magnitude  between  algorithms  is  calculated.  By  noting  ANIK-Fl’s 
intensity  at  Pd  -  0.5  for  each  algorithm,  the  apparent  magnitude,  AM,  between  them  can 
be  determined. 


AM  =  -2.51og10(^— (3.30) 

\I  ref  —  B) 

I  is  the  total  intensity  in  the  pixel  containing  ANIK-F1  when  either  the  equal  or 
unequal  cost  algorithm  crosses  the  50  percent  detection  threshold,  and  Jre f  is  the  total 
intensity  when  the  current  BHT  crosses  50  percent.  A  larger  positive  AM  indicates  a 
dimmer  object  in  comparison  to  the  reference  object. 


Table  3.3:  Apparent  magnitude  improvement  of  Equal  and  FTnequal  MHT  over  the  current 
BHT  at  50  percent  detection  threshold. 


13  March 

14  March 

15  March 

ECEP  MHT 

1.67 

0.46 

1.06 

UCEP  MHT 

2.40 

0.87 

1.44 
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On  all  three  nights  both  of  the  MHT  algorithms  used  saw  an  increase  in  AM  over  the 
BHT.  In  addition,  there  is  an  improvement  of  0.73,  0.41,  and  0.38  in  apparent  magnitude 
when  the  UCEP  algorithm  is  used. 

The  detection  performance  at  the  current  SST  program  specified  Pf  is  important,  but 
does  not  completely  compare  the  detection  performance  against  a  range  of  false  alarm  rates. 
To  investigate  Pd  performance  across  varying  levels  of  Pf ,  ROC  curves  are  generated  at 
multiple  intensity  intensity  levels. 

3.4.2  Receiver  Operating  Characteristic. 

A  ROC  curve  is  a  method  to  demonstrate  Pd  and  Pf  as  a  function  of  the  selected 
threshold  levels.  To  generate  a  ROC  curve,  sets  of  Pd  and  Pf  data  pairs  are  calculated 
for  a  specific  threshold  level  r.  The  threshold  is  varied  until  a  full  range  of  detection  and 
false  alarm  probabilities  between  0  and  1  are  calculated.  The  ROC  curve  for  the  equal  cost 
method  is  calculated  with  the  following  relation. 


Pe(j)  =  ^  FMp  +  T)Pr,(Vp)&tip  (3-31) 

vP 

The  threshold  can  be  combined  with  rjp  in  various  ways,  but  addition  is  chosen  because 
it  allows  for  an  easy  way  to  determine  r  values  needed,  and  it  works  even  if  the  value  for 
T]p  is  zero.  Values  for  the  threshold  r  are  chosen  to  ensure  two  key  factors.  The  first  is  that 
all  possible  values  of  Pf  and  Pd  are  obtained.  The  second  is  that  there  are  enough  points  on 
the  sampling  grid  to  generate  a  smooth  curve.  The  unequal  cost  ROC  curve  is  found  with 
the  following  equation. 

r„+r 

Pu(j)=  J  p(Lu)dLu  =  F£ii(Tu  +  T)  (3.32) 

— oo 

Due  to  the  large  quantity  of  pixels  being  tested,  a  low  Pf  is  required  to  give  a 
manageable  number  of  false  alarms  across  each  image.  Due  to  the  importance  of  low 
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Pf  values,  the  ROC  curve  is  plotted  as  a  semi-log  plot.  Taking  the  log  of  Pf  provides  better 
insight  into  algorithm  performance  across  a  larger  set  of  values. 

Figure  3.5  includes  ROC  curves  for  both  the  ECEP  and  UCEP  MHT  algorithms.  The 
ROC  curves  are  produced  with  15  frame  windows  centered  at  three  different  frames  from 
the  data  collection  on  13  Mar  2012. 
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Figure  3.5:  ROC  curves  centered  at  frames  (A)  800  (B)  850  (C)  900  from  13  March  2012 
using  a  15  frame  window  with  M  -  10. 


In  (A),  the  window  is  centered  at  frame  800.  At  this  point  in  the  collection,  the 
satellite  is  still  very  visible.  The  ECEP  and  UCEP  algorithms  both  achieve  a  greater  than 
80  percent  Pd  at  up  to  10e-10  Pf.  The  ECEP  algorithm  drops  off  more  quickly  than  the 
UCEP  algorithm,  giving  up  to  an  80  percent  increase  when  Pf  =  lOe  -  25.  In  (B),  frame 
850  is  used  as  the  center  frame.  The  satellite  is  more  difficult  to  detect  in  this  window. 
The  UCEP  algorithm  performs  much  better  for  a  large  variety  of  Pf  values,  including  up 
to  a  100  percent  difference  at  10e-12.  In  (C),  a  window  around  frame  900  is  utilized.  The 
satellite  is  very  difficult  to  detect  at  this  low  intensity.  Looking  at  an  image  at  this  point,  it 
would  be  hard  to  determine  if  an  object  is  present  or  not,  as  shown  in  Figure  3.6.  Figure 
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3.6  shows  ANIK-F1  at  three  points  as  it  enters  eclipse  on  13  March.  Both  algorithms  have 
a  lower  probability  of  detection  for  the  same  probability  of  false  alarm  when  compared 
to  earlier  frames.  At  the  realistic  Pf  values  that  might  be  used  in  the  SST  system,  both 
algorithms  are  below  10  percent  detection. 


(A)  (B)  (C) 


Figure  3.6:  ANIK-F1  at  3  points  during  the  eclipse,  captured  with  the  SST.  The  frames 
used  are  (A)  800,  (B)  850,  and  (C)  900  on  13  March.  ANIK-F1  loses  intensity  as  it  enters 
the  eclipse  and  more  close  resembles  the  background. 


The  next  factor  investigated  is  the  number  of  hypotheses,  M,  used  in  the  test.  As 
described  in  section  3.2.1,  the  number  of  hypotheses  used  may  impact  the  detection 
performance  of  the  algorithm.  To  compare  M,  the  UCEP  structure,  £„,  is  used.  The 
data  collected  on  13  Mar  2012  is  analyzed  in  this  chapter  with  other  nights  demonstrating 
similar  performance.  Figure  3.7  shows  ROC  curves  for  M  values  of  2,  6,  and  10,  for 
varying  satellite  intensities.  For  M  =  6  four  comer  sub-pixel  locations  are  used  along  with 
the  center  of  the  pixel.  To  generate  these  corner  locations,  both  a  and  a>  shifts  of  ±30//m 
are  used.  For  M  -  10,  top,  bottom,  and  side  locations  are  added.  These  locations  use  a 
shift  of  ±30// m  in  either  a  or  to,  and  no  shift  in  the  other  dimension. 
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Figure  3.7:  ROC  curve  using  UCEP  MHT,  X„,  at  frame  (A)  800  (B)  850  and  (C)  900  from 
13  March  2012  using  a  15  frame  window.  Three  values  of  M  are  compared:  2,  6,  and  10. 


Several  interesting  trends  are  apparent  as  the  value  of  M  changes.  One  important 
observation  is  that  the  number  of  hypotheses  used  has  an  impact  on  the  performance  of  the 
algorithm.  This  effect  is  more  noticeable  for  higher  intensity  objects.  In  (A),  the  satellite 
is  still  intense,  and  the  M  -  2  ROC  curve  is  lower  than  both  the  6  and  10  hypotheses 
case.  The  performance  difference  is  an  up  to  60  percent  lower  detection  rate  for  the  same 
false  alarm  rate,  implying  that  the  MHT  performs  better  than  the  BHT.  As  the  object 
loses  intensity  and  enters  the  eclipse,  the  separation  between  the  MHT  and  the  BHT  is 
decreased.  In  (B),  at  frame  850  the  max  detection  separation  is  approximately  20  percent, 
and  in  (C)  at  frame  900  the  separation  is  approximately  15  percent.  Another  interesting 
result  is  that  the  difference  between  the  number  of  hypotheses  beyond  the  binary  case  is 
less  significant.  Using  M  -  6,  the  detection  performance  is  never  below  the  M  -  10  case. 
The  improvement  ranges  from  0-5  percent  depending  on  the  night  being  processed  and  the 
false  alarm  rate.  This  result  demonstrates  that  there  is  a  trade-off  between  increasing  the 
number  of  sub-pixel  locations  to  capture  potential  aliased  PSF  and  creating  additional  false 
alarms.  The  choice  of  M  =  6  is  due  to  the  fact  that  it  is  the  first  multiple  hypothesis  case 
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that  will  symmetrically  cover  potential  sub-pixel  locations.  This  is  done  by  choosing  the 
four  corners  and  the  center. 

One  important  thing  to  note  is  that  the  M  =  2  case  shown  in  Figure  3.7  is  a  BHT,  but  is 
not  the  point  detector  currently  used  by  SST.  There  are  two  differences  between  this  binary 
correlation  and  a  point  detector.  The  first  difference  is  that  in  the  correlation,  a  15x15  pixel 
window  is  used  to  make  the  decision  between  H0  and  Hi,  and  only  a  single  pixel  is  used  in 
the  point  detector  .  Secondly,  the  correlation  uses  outlier  removal  as  described  in  section 
3.4. 

In  addition  to  the  performance  increase,  using  M  =  6  also  reduces  the  required 
computations  to  calculate  £„.  By  using  six  hypotheses,  only  five  alternate  hypothesis 
terms  need  to  be  calculated,  compared  to  nine  when  M  =  10,  reducing  computation  cost 
by  approximately  45  percent,  assuming  the  processing  times  for  each  hypothesis  are  equal. 
The  selection  of  M  gives  a  way  to  control  the  computational  time  required  for  the  unequal 
cost  method,  but  in  general  this  method  is  more  costly  than  the  currently  implemented 
point  detector.  Due  to  the  fact  that  £u  is  not  linear,  processing  techniques  like  Fourier 
transforms  cannot  be  used  to  increase  processing  speed.  Currently  the  point  detector  in  the 
SST  system  runs  in  real  time,  and  can  create  detections  as  quickly  as  images  are  captured. 
Using  the  new  proposed  UCEP  MHT,  the  same  processing  technique  will  not  be  possible. 
This  allows  the  user  a  choose  between  processing  speed  and  detection  performance. 

3.5  Full  Frame  Implementation 

The  UCEP  algorithm  presented  in  this  research  is  more  computationally  complex  than 
previously  proposed  methods.  To  implement  the  algorithms  into  the  current  SST  data 
processing  pipeline,  it  is  important  to  investigate  the  feasibility  of  achieving  real  time  or 
near  real  time  analysis.  Looking  at  the  final  detection  algorithm  for  the  unequal  cost  method 
gives  insight  into  possible  methods  for  processing  the  data  more  quickly. 
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The  most  computationally  expensive  portion  of  the  algorithm  is  performing  operations 
on  the  large  matrices.  As  mentioned  previously,  the  full  frame  data,  2x2  binned,  contains 
6144x4096  pixels.  The  first  step  to  increase  the  processing  speed  is  to  separate  the 
argument  of  the  exponential  function  in  X„  into  data  dependent  and  non-data  dependent 
calculations. 


jgQ  "y  QZ.  ^ x  *'y 

hk(x, y)D(x, y) - T  E  E  hk(yX ’ y ^  ~  E  E  hk(yX ’ (3-34) 

x=\  y=  1  X—  1  >•=!  X—  1  y=l 

^  .  -Nx"  . 
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ip  represents  the  non-data  dependent  portion  of  the  equation.  Both  B  and  cr  are 
background  data  statistics  that  need  to  be  computed  in  each  frame,  and  are  needed  for 
both  ECEP  and  UCEP  algorithms.  The  only  variable  changing  from  frame  to  frame  is  the 
data,  D(x,y).  The  significant  difference  is  the  correlation  between  the  data  and  the  kth  PSF. 
The  UCEP  algorithm  has  two  additional  data  processing  steps;  the  exponential  function 
and  a  sum  of  the  M  —  1  total  hypotheses. 

Allowing  the  window  of  Nx  and  Ny  to  go  to  the  entire  frame  allows  for  simultaneous 
processing  of  an  entire  data  frame.  A  small  window  no  longer  needs  to  be  used  to  cycle 
through  the  image,  decreasing  the  processing  time.  MATLAB  has  a  function  that  allows 
data  processing  on  the  graphics  card.  Using  this  function,  the  correlation  can  be  passed  to 
a  compatible  graphics  card.  The  correlation  function  can  be  computed  more  quickly  with 
this  method,  allowing  for  a  faster  overall  processing  time. 

The  full  frame  data  is  processed  with  a  HP  Z820  workstation.  The  processor  in  the 
workstation  is  a  Intel  Xeon  CPU  E5-2650  at  2.6GHz  with  64  GB  of  RAM.  The  software 
used  is  Windows  7  and  MATLAB  R2015b.  The  graphics  processor  used  is  a  NVIDIA 
Quadro  K6000.  Defining  TEcup  and  TECep  as  the  time  it  takes  the  computer  described  to 
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process  the  respective  algorithms,  the  average  and  maximum  processing  times  are  shown 
below. 


Table  3.4:  Processing  times  for  equal-  and  unequal-cost  algorithms  processed  on  graphics 
card. 


T ECEP 

T ECUP 

AT 

Mean 

187  ms 

237  ms 

26.7% 

Maximum 

408  ms 

568  ms 

39.2% 

AT  represents  the  additional  processing  time  required  by  the  ECUP  algorithm  as 
a  percentage.  As  Table  3.4  demonstrates,  there  is  a  definite  jump  in  processing  time 
required  by  the  UCEP  algorithm.  The  26.7  percent  additional  time  may  or  may  not  impact 
real  time  operation  of  the  processing  depending  on  the  required  time  for  other  processes. 
These  processes  include  loading  the  data  frames,  clustering,  orbit  determination,  and  other 
post  detection  processing.  These  processes  are  common  to  both  the  ECEP  and  UCEP 
algorithms. 

3.6  Conclusions 

This  chapter  presents  a  method  of  increasing  the  detection  capability  of  the  SST  over 
the  current  BHT  algorithm.  Due  to  the  spatial  aliasing  caused  by  the  pixel  size,  a  MHT 
is  used  to  account  for  potential  different  sub-pixel  shifts  resulting  in  changes  to  the  shape 
and  distribution  of  the  PSF.  These  PSFs  are  used  to  correlate  with  received  data  to  detect 
objects.  It  is  demonstrated  that  the  probability  of  detection  is  increased  by  using  an  UCEP 
MHT  detection  algorithm  when  compared  to  an  ECEP  algorithm.  Pd  gains  of  up  to  80 
percent  are  observed  for  the  same  P  f  rate.  The  number  of  hypotheses  used  within  the  MHT 
is  also  investigated.  It  is  found  that  there  is  a  relationship  between  the  number  of  hypotheses 
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used  and  detection  performance.  The  results  show  that  a  MHT  clearly  outperforms  a  BHT, 
and  that  using  6  hypotheses  slightly  outperformed  10,  with  the  additional  benefit  of  reduced 
computational  time.  With  the  M=6  hypothesis  assumption,  the  aliased  PSFs  proposed  more 
accurately  match  the  observed  PSF  in  the  data  without  raising  the  probability  of  false  alarm 
by  introducing  additional  hypotheses. 

The  key  contributions  of  this  research  are  a  new  reduced  algorithm  for  implementing 
an  UCEP  MHT  based  on  received  data  with  a  Gaussian  distribution,  and  a  method  for 
determining  average  Pd  and  Pj  based  on  a  general  cost  structure.  This  chapter  investigates 
the  limiting  cases  of  cost  selections.  The  first  being  that  all  errors  have  the  same  cost, 
resulting  in  a  large  emphasis  being  placed  on  correctly  identifying  the  correct  sub-pixel 
location.  The  second  limiting  case  is  the  UCEP  algorithm.  The  UCEPdoes  not  put  any 
cost  on  errors  in  sub-pixel  location,  and  instead  uses  that  data  in  essence  as  a  nuisance 
parameter.  The  ability  to  generate  ROC  curves  for  any  intermediate  case  can  be  useful  in 
determining  the  best  cost  structure  for  a  specific  application. 

There  are  several  areas  of  future  research  that  could  improve  or  add  to  the  results 
presented.  This  research  focuses  on  received  data  that  is  assumed  to  be  Gaussian. 
Through  this  research,  a  method  for  implementing  a  different  model  for  received  data  is 
developed  that  can  be  used  on  any  potential  distribution  of  received  data.  The  limiting 
cases  of  costs  are  considered  in  this  chapter.  There  are  many  additional  cost  structures 
in  between  the  two  presented  here.  For  specific  applications  these  cases  may  give  better 
detection  performance.  The  drawback  is  that  they  will  not  reduce  as  well  as  the  equal 
or  unequal  cost  schemes  presented  here.  Additionally,  the  unequal  cost  method  here  is 
more  computationally  complex  than  the  currently  implemented  algorithms.  Simplification 
and  efficient  implementation  are  not  a  focus  of  this  research.  To  use  this  algorithm 
operationally,  these  factors  will  need  to  be  investigated  further. 
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In  the  next  chapter  a  Bayes  risk  MHT  is  further  investigated  by  developing  a  detection 
algorithm  with  unequal  prior  probability. 
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IV.  Investigating  Multiple  Hypothesis  Test  Prior  Probabilities 


4.1  Introduction 

The  main  goal  in  this  chapter  of  the  dissertation  is  to  create  a  detection  algorithm 
that  uses  an  unequal  prior  probability  assignment  to  improve  detection  capabilities.  This 
chapter  addresses  the  research  question:  Do  the  assignments  of  a  priori  probabilities  in  a 
MHT  improve  the  detection  performance? 

The  distribution  of  space  objects  on  the  CCD  of  a  ground-based  telescope,  and  how 
that  sub-pixel  position  translates  to  several  potential  aliased  PSFs  is  key  in  answering  this 
research  question.  Due  to  the  large  number  of  potential  space  objects,  their  arbitrary 
location  relative  to  the  observation  point,  and  physical  properties  of  the  optics,  it  is 
acceptable  to  assume  that  a  RSO  has  a  uniform  probability  of  occurring  anywhere  within 
a  single  pixel.  The  question  to  be  answered  in  this  chapter  is:  assuming  a  uniform 
distribution  across  a  single  CCD  pixel,  what  are  the  resulting  prior  probabilities,  nk,  in 
a  MHT?  Additionally,  if  a  detection  algorithm  is  derived  using  these  probabilities,  is  it  able 
to  outperform  previously  proposed  hypothesis  tests? 

As  described  in  Chapter  2,  due  to  aliasing  present  in  many  telescopes,  the  sub-pixel 
position  can  affect  both  the  shape  and  location  of  the  resulting  PSF.  This  chapter  covers 
several  different  topics  relating  to  the  investigation  of  the  prior  probabilities.  These  aspects 
include:  background  and  motivation  for  this  type  of  MHT,  a  decision  space  analysis  for 
assigning  prior  probabilities,  optical  modeling  theory,  and  the  development  of  the  detection 
algorithm.  The  algorithm  is  tested  against  both  simulated  space  objects  and  data  collected 
from  the  S ST. 
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4.2  Theory 

To  implement  an  unequal-prior  detection  algorithm,  there  are  two  important  factors. 
One  factor  is  the  ability  of  the  algorithm  to  reduce  to  a  calculable  statistic.  Also  pertinent 
is  how  well  the  algorithm  detections  objects,  and  realistic  and  representative  values  for  the 
prior  probabilities  nk.  This  section  covers  the  theory  on  both  of  these  topics. 

4.2.1  Optical  Model. 

To  mitigate  the  effects  of  aliasing,  candidate  PSFs  are  needed  for  the  MHT  to  compare 
against.  There  are  two  approaches  to  producing  these.  The  first  is  to  capture  representative 
PSFs  from  the  collected  data.  The  drawback  to  this  method  is  that  the  PSF  is  specific  to 
the  conditions  of  the  night  captured,  including  the  seeing,  rQ,  and  the  sub-pixel  location  of 
the  object.  A  large  sampling  of  objects  would  potentially  need  to  be  investigated  in  order 
to  obtain  complete  models  for  the  M  alternate  hypotheses. 

Alternatively,  this  research  uses  an  optical  model  to  create  representative  PSF  based 
on  the  input  conditions.  A  more  detailed  derivation  is  described  in  section  3.1.  The 
final  PSFs  hk(x,y),  shown  in  equation  (3.10),  are  created  through  a  combination  of  four 
important  factors:  the  optical  PSF  including  lens  aberrations  and  telescope  parameters,  a 
long  exposure  atmosphere  model  [17],  any  shift  from  the  center  of  the  CCD,  and  the  spatial 
aliasing  effect  modeled  as  a  blurring  function. 

The  modeled  PSFs  hk(x,y )  are  used  for  several  purposes  in  this  chapter.  It  is  used  to 
evaluate  the  decision  space  and  determine  prior  probabilities,  for  data  creation,  and  is  a  key 
factor  in  the  detection  algorithm.  Next,  the  decision  space  for  the  hypotheses  in  this  test  is 
investigated  to  determine  the  prior  probability. 

4.2.2  Bayes  Cost  and  Priors  Discussion. 

To  make  a  decision  if  an  object  is  present,  a  Bayes  Criterion  is  used  [25].  This  equation 
provides  a  method  of  choosing  the  hypothesis  that  results  in  less  risk  K  on  average. 
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(4.1) 


M- 1  M- 1  r 

^Cifc  I  p(D|ff*)dD 

i=0  £=0  ^Z' 

M  is  the  total  number  of  hypotheses  considered.  The  cost,  Clk ,  represents  the  impact  of 
choosing  hypothesis  i  when  hypothesis  k  has  occurred.  The  value  of  Clk  ranges  from  zero 
to  one,  with  one  giving  the  largest  cost  penalty  to  a  decision.  The  prior  probabilities,  nk, 
represent  how  likely  each  potential  hypothesis  is  to  occur.  The  priors  range  between  zero 
and  one,  and  must  sum  to  a  total  of  one.  When  using  a  Bayes  Criterion-based  MHT,  there 
are  several  variables  that  impact  the  algorithm  that  results  from  the  hypothesis  test.  The 
two  discussed  in  this  paper  are  the  costs,  Cik  and  the  a  priori  probabilities,  nk.  D  is  a  matrix 
containing  one  single  frame  of  data.  Hk  is  the  kxh  hypothesis.  H{)  is  the  null  hypothesis, 
where  it  is  assumed  that  no  space  object  is  present.  Hi  through  HM_ i  are  the  alternate 
hypotheses.  These  correspond  to  instances  where  the  space  object  is  considered  present, 
and  each  hypothesis  signifies  a  different  sub-pixel  position.  These  sub-pixel  positions  are 
described  in  Table  4.1. 

Different  cost  and  a  priori  probability  approaches  have  been  researched,  and  there  are 
several  differences  between  the  approaches.  The  first  area  of  difference  is  the  assignment 
of  cost.  An  Equal-Cost  Equal-Prior  (ECEP)  test,  originally  proposed  in  [24],  penalizes  the 
algorithm  for  incorrectly  deciding  between  two  alternate  hypotheses.  This  emphasis  on 
determining  the  correct  sub-pixel  position  can  lead  to  more  accurate  sub-pixel  position 
estimates,  potentially  at  the  cost  of  the  detection  of  space  objects.  Alternatively,  an 
Unequal-Cost  Equal-Prior  (UCEP)  assumption  will  not  penalize  for  selecting  the  wrong 
alternate  hypothesis,  as  long  as  an  object  is  correctly  detected  [67].  Unequal-cost 
algorithms  are  more  computationally  complex  and  do  not  reduce  to  a  SNR  sufficient 
statistic.  Unequal-cost  algorithms  are  not  discussed  further  in  this  chapter. 

The  other  variables  where  different  assumptions  can  be  made  in  a  Bayes  Risk  MHT 
are  the  a  priori  probabilities.  In  this  chapter,  two  different  a  priori  probability  approaches 
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are  investigated.  They  are  an  Equal-Cost  Equal-Prior  (ECEP)  and  an  Equal-Cost  Unequal- 
Prior  (ECUP).  ECEP  assumes  the  following  cost  and  probability  assignment: 


f  0,  i  =  k 
I  1,  i  4  k 


(4.2) 


The  ECUP  algorithm  assumes  the  same  cost  assignment,  but  the  prior  probabilities, 
TTjfc,  for  each  hypothesis  k  are  not  necessarily  the  same  or  j~r  The  prior  probabilities  are 
typically  assumed  to  be  equal  if  there  is  no  prior  knowledge  of  the  system  or  method  of 
determining  accurate  values  for  This  chapter  considers  several  methods  for  finding 
accurate  prior  probability  values,  and  determines  if  more  binary  detections  can  be  found 
with  the  resulting  algorithm. 

4.2.3  Decision  Space  Analysis. 

A  key  element  to  creating  an  algorithm  that  takes  into  account  an  unequal  prior 
probability  assumption  is  determining  accurate  prior  probabilities  for  each  sub-pixel 
hypothesis,  Hk.  Without  a  method  of  accurately  determining  these  values,  the  actual 
benefit  of  an  ECUP  test  may  not  be  realized.  The  assignment  of  priors  reflects  how  point 
sources  present  in  different  locations  within  a  CCD  translate  into  potential  PSFs.  This  paper 
investigates  three  potential  methods  for  segmenting  the  pixel  into  a  decision  space.  These 
methods  are:  a  distance-based  metric,  a  correlation  metric,  and  an  empirical  method.  Each 
method  is  described  and  analyzed  in  the  following  paragraphs. 

A  critical  assumption  common  with  all  of  these  methods  is  that  a  space  object  is 
equally  likely  at  any  position  within  a  pixel.  This  sub-pixel  position,  (a,  to),  is  the  physical 
distance  in  pm  within  a  2x2  binned  pixel.  This  assumption  implies  a  uniform  probability 
across  the  entire  pixel.  To  segment  the  decision  space,  a  sub-pixel  map  is  created.  In  the 
map,  each  position  tested  has  the  same  probability  of  occurrence.  To  ensure  a  fine  enough 
coverage  of  the  entire  pixel,  sub-pixel  positions  are  tested  every  \pm  in  both  a  and  to.  This 
results  in  31  positions  in  each  dimension  and  961  positions  total. 
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In  this  chapter,  M  =  10  hypotheses  are  used  in  the  test.  These  hypotheses  are  PSFs 
generated  through  the  optical  model  where  a  space  object  is  assigned  to  be  in  a  defined 
location  within  a  pixel.  As  previously  mentioned,  the  hypotheses  positions  are  defined  by 
a  and  to,  the  position  within  a  30x30/vm  pixel,  with  the  origin  being  the  center  of  the  pixel. 
The  locations  within  a  pixel  are  the  corners  where  a  =  ±15  and  to  =  +15,  sides  where 
a  =  0  and  co  =  ±15,  top  and  bottom  where  a  =  +15  and  to  =  0,  and  the  center  of  each 
pixel. 


Figure  4.1:  Decision  space  within  a  single  CCD  showing  all  sub-pixel  locations,  along  with 
the  locations  of  the  nine  alternate  hypotheses. 


There  are  multiple  reasons  for  using  these  positions.  First,  they  give  an  even 
representation  of  the  potential  resulting  PSFs  based  on  the  coverage  of  the  entire  pixel. 
Additionally,  they  can  reduce  computational  complexity  and  reduce  the  number  of  tests, 
due  to  sharing  the  corners  and  sides  hypotheses  between  adjacent  pixels.  In  this  case, 
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rather  than  testing  all  nine  alternate  hypotheses,  only  five  are  needed  for  each  pixel.  This 
approach  was  first  proposed  in  [24].  An  alternate  layout  of  hypothesis  positions  with  M  =  6 
is  proposed  in  chapter  3,  but  is  not  investigated  in  this  chapter.  M  -  10  hypotheses  are  used 
in  this  research  to  compare  against  the  previously  proposed  M  =  10  ECEP  algorithm. 

The  first  method  considered  is  the  correlation  metric.  Both  the  ECEP  and  ECUP 
algorithms  are  based  on  a  matched  filter,  or  correlation  test.  By  determining  which 
hypothesis  correlates  most  closely  with  each  sub-pixel  position  tested,  an  assignment 
matrix  can  be  formed  by  noting  the  most  closely  correlated  hypothesis  at  each  sub-pixel 
position  a,  co.  This  is  done  with  the  following  equation: 


'Hc(o',  a>)  =  argmax 


NX  Ny 

XX 

x=l  y=l 


Ta,oj(x,y)hk(x,y ) 
(Tk 


(4.3) 


crk  is  a  normalization  term  based  on  the  kth  hypothesis  present  in  the  ECEP  algorithm, 

/v,  Ny 

crk  -  J  X  X  K(x,y).  (Hc(a,co)  is  an  entry  in  the  hypothesis  assignment  matrix  T£r. 

\  v  i  y  : 

Each  coordinate  represents  a  sub-pixel  position  and  contains  the  closest  hypothesis  match. 
Ta,o(x,y)  is  a  16x16  pixel  modeled  PSF  based  on  a  space  object  at  the  sub-pixel  location 
(a,  cl>).  An  important  note  is  that  the  PSFs  for  T,rM(x,y)  and  hk(x,y )  are  generated  by 
performing  the  appropriate  shift  on  the  highly  sampled  model  PSFs  and  downsampled  as 
described  in  equation  (3.9). 

The  second  proposed  method  for  assigning  a  sub-pixel  position  to  a  hypothesis  being 
considered  is  a  distance-based  metric.  The  distance  vector  in  this  case  is  between  the  kth 
hypothesis  position  and  the  proposed  a,  ft  sub-pixel  position.  There  are  several  methods 
for  measuring  the  distance,  or  size  of  a  vector.  In  this  research,  a  2-norm  is  used  to  find 
the  “closest  point”  between  the  sub-pixel  position  being  tested  and  the  locations  of  the 
hypotheses. 


<Hd(a,  to)  =  argmin(||i4||) 


(4.4) 
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Vk  is  a  vector  between  each  hypothesis  k  located  at  (ak,  tok),  and  the  sub-pixel  position, 
(or,  to),  being  tested.  'HJ.a,  to)  is  the  hypothesis  assignment  matrix  for  the  distance-based 
metric.  At  each  sub-pixel  position,  nine  vectors  are  created  and  the  smallest  2-norm  is 
selected  as  the  hypothesis  that  best  represents  the  sub-pixel  position  being  tested. 

Figure  4.2  shows  the  decision  space  analysis  for  a  single  CCD  pixel  resulting  from 
the  correlation  and  distance-based  metrics.  The  plots  demonstrate  the  physical  layout  of 
the  decision  space  and  the  boundaries  between  each  hypothesis.  Each  sub-pixel  position  is 
assigned  a  corresponding  hypothesis  and  is  grouped  into  a  section  with  similarly  assigned 
hypotheses. 


to  Position  (p.m)  to  Position  (fim) 


Figure  4.2:  Decision  space  analysis  for  a  single  SST  CCD  pixel.  The  hypothesis 
assignment  matrix  for  (A)  the  correlation  metric,  'Hc,  and  (B)  distance  metric, Each 
color  shade  corresponds  to  a  similar  hypothesis.  Sub-pixel  positions  are  tested  every  1/vm. 


Looking  at  the  distribution  of  hypothesis  assignments  in  Figure  4.2,  there  are 
similarities  between  the  correlation  and  distance-based  metrics.  Both  segment  the  decision 
space  into  rectangles  based  on  the  location  of  the  hypotheses.  The  correlation  metric  has 
sub-pixel  positions  between  H5  and  //x,  as  well  as  H4  and  H1,  that  do  not  directly  create 
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perfect  rectangular  boundaries.  These  are  due  to  optical  effects  not  included  in  the  distance 
metric,  such  as  lens  or  mirror  aberrations  and  atmospheric  effects. 

The  final  method  being  considered  to  determine  prior  probabilities  is  an  empirical 
method.  The  empirical  method  uses  detection  made  on  collected  SST  data  to  determine 
how  often  each  hypothesis  is  observed.  This  method  depends  on  implementing  the  ECEP 
algorithm  and  noting  the  selected  hypothesis.  One  drawback  is  that  the  empirical  method 
only  analyzes  objects  that  it  detects.  The  inherent  assumption  when  using  this  method  is 
that  non-detectable  objects  will  have  the  same  spatial  distribution  as  detectable  objects. 
Using  the  ECEP  algorithm,  the  following  SNR  equation  can  be  used  to  assign  sub-pixel 
positions  [24]: 


Nx  Ny 


SNR*  =  Z  £ 


(DXo,yo(x,y)  -  B)  hk(x,y) 


X=\  V=1 


CT 


crk 


(4.5) 


DXo,y0  (x,y)  is  a  data  window  that  is  Nx  by  Ny  pixels  centered  at  x0,y0. 


'HeiXo,  y0)  =  argmax  [SNRk]  (4.6) 

k 

The  three  methods  described  to  this  point  all  have  an  associated  hypothesis  assignment 
matrix,  'H.  This  matrix  gives  each  sub-pixel  position  the  hypothesis  that  it  most  closely 
resembles.  The  prior  probability  values  calculated  from  all  three  methods  are  included  in 
section  4.4.  One  distinction  between  the  methods  is  that  'Hc  and  'H(/  represent  a  single  CCD 
pixel,  while  HLe  represents  an  entire  frame  of  SST  data.  Therefore,  H.,  does  not  provide 
insight  into  how  the  decision  space  inside  the  pixel  is  physically  distributed  as  shown  in 
Figure  4.2,  only  the  prior  probability  of  each  hypothesis  occurring. 

4.2.4  Hypothesis  Test  Derivation. 

The  next  step  is  to  implement  the  potential  nk  values  discussed  in  the  previous 
section  into  a  detection  algorithm  and  determine  how  the  change  in  assumptions  affects 
the  detection  performance  of  the  algorithm.  As  mentioned  previously,  with  an  unequal- 
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prior  approach  there  are  two  choices  for  costs:  Equal-Cost  Unequal-Prior  (ECUP)  and 
Unequal-Cost  Unequal- Prior  (UCUP).  When  implementing  an  equal-cost  assumption,  the 
Bayes  criterion  [25]  can  be  reduced  to  the  following  relation: 


max 

k 


7Tk  p(D\Hk) 


%  1 


(4.7) 


n0  p(D\H0 ) 

where  p(D\Hk)  is  the  conditional  PDF  of  the  data,  given  that  hypothesis  k  is  true.  D  is  a 
matrix  that  represents  a  received  image.  The  hypothesis  test  will  select  the  max  ratio  of 
probabilities  and  conditional  PDFs  and  compare  the  resulting  ratio  against  one.  Based  on 
the  result  of  the  comparison,  a  decision  on  which  hypothesis  is  selected  can  be  made.  If  the 
ratio  is  less  than  one,  it  is  assumed  that  no  object  is  present.  If  the  ratio  is  greater  than  one, 
an  object  is  assumed  present  and  the  largest  hypothesis  k  is  selected.  A  varying  selection 
of  costs  will  result  in  a  different  ratio  test.  Alternate  cost  assumptions  are  not  investigated 
in  this  chapter. 

To  reduce  equation  (4.7),  a  model  for  what  the  received  data  D  contains  is  needed. 
This  model  is  described  in  equation  (4.8). 


f  0hk(x,y)  +  B  Hh,k  6  [1,9] 

I(x,y)  =  (4.8) 

[  B  Hq 

6  is  the  modeled  RSO  intensity.  hk(x,y)  is  the  PSF  based  on  the  kth  hypothesis 
described  in  section  4.2.1  and  B  is  the  background  photo  count.  I(x,y )  represents  the 
average  expected  value,  but  there  will  be  noise  present  in  the  data.  As  a  result,  a  noise 
model  is  needed  to  completely  represent  the  conditional  PDF.  In  this  case,  the  model 
assumes  Poisson  counting  as  the  dominant  noise  source.  In  previous  research  efforts,  both 
Poisson  and  Gaussian  noise  assumptions  have  been  made.  One  reason  for  this  is  that  it 
has  been  shown  that  for  a  flat  background,  as  is  the  case  in  this  situation,  both  Gaussian 
and  Poisson  assumptions  lead  to  the  same  algorithm  [22,  35].  Multiple  sources  of  noise 
tend  to  give  an  overall  Gaussian  distribution  to  the  noise  present  in  the  system.  Alternative 
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noise  assumptions  are  investigated  further  in  Chapter  5  [68].  Assuming  an  independent 
joint  Gaussian  noise  distribution,  the  PDF  ratio  in  equation  (4.7)  results  in  the  following 
relation: 


NxNy 


(  NX  Ny 


p(T>\Hk)  =  (2 ncr2}  2  exp  -  ^  ^ 

V  -Y—  1  V=  I 


(D(x,y)  -  I(x,y))2 
2  cr2 


(4.9) 


Nx  and  Ny  are  the  number  of  pixels  in  the  x  and  y  direction,  cr  is  the  measured  standard 
deviation  of  the  received  data  D(x,y).  Combining  equations  (4.7)  and  (4.8),  combining 
exponentials,  and  expanding  gives  the  following  equation: 


exp 


NX  Ny 
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Hq  7tk 
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This  equation  can  be  further  reduced  by  taking  the  natural  log,  removing  the 
exponential  function. 


NX  Ny 


(02/^2(- x>y)  -  2 6h(x,y)(D(x,y)  -  B ))  ^  In 
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x=l  y=l 

It  is  obvious  that  for  an  equal-prior  assumption,  the  natural  log  of  prior  probabilities 
goes  to  zero.  This  assumption  gives  one  less  term  on  the  threshold  side  of  the  equation  to 
be  concerned  with.  Continuing  to  reduce  the  left  side  of  the  equation  gives: 
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The  left  side  of  equation  (4.12)  now  resembles  a  SNR,  but  without  the  standard 
deviation  of  the  data.  To  give  a  true  SNR,  both  sides  are  divided  by  cr  to  give  the  following: 


|ln 

,  ,  ^  Ho 

X=l  y= 1  u 


Nx  Ny 

6  I  I  h2(x,  y) 

cr  x=\  >=i 
6  2  cr 


(4.13) 


The  left  side  of  equation  (4.13)  contains  a  correlation  of  the  received  data  D(x,y )  and 
the  kth  PSF.  This  implies  that  the  algorithm  is  looking  for  how  similar  the  received  data  and 
the  expected  PSF  are.  Introducing  a  normalization  factor  to  ensure  that  the  SNR  statistic  is 
unit  variance  [24]  gives  the  following  equation: 


NX  Ny 
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(4.14) 


wk  r 

where  crk  is  a  normalization  factor  to  keep  the  SNR  a  zero  mean  unit  variance  random 
variable  [24].  Investigating  the  threshold  gives  insight  into  the  effect  on  the  threshold  of 
unequal  prior  probability. 

Combining  equations  (4.7)  and  (4.20)  and  grouping  terms  dependent  on  K  gives  this 
reduced  equation: 


max  [SNR/;  -  Wk\  =  SNRWk  %  Yk.  (4.15) 

k  H0 

where  Wk  is  a  weighting  factor  based  on  the  non-equal  prior  probabilities.  If  n0  =  nk,  then 
the  weighting  term  will  reduce  to  zero  and  equal  the  threshold  in  the  equal  prior  probability 
case.  The  ratio  of  7r0  to  nk  can  increase  or  decrease  the  weighting  factor,  changing  the 
probability  of  false  alarm  and  detection.  Increasing  the  probability  of  any  hypothesis  other 
than  the  null  hypothesis  will  lower  the  effective  threshold  value. 
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Modifying  the  SNR  by  the  weighting  factor  changes  the  mean,  and  therefore  the 
algorithm  does  not  have  exactly  the  same  Pj  as  [24].  In  this  research  the  Pf  is  determined 
through  applying  the  algorithm  to  simulated  data  containing  no  RSO.  Due  to  this  factor, 
a  false  alarm  analysis  is  also  included,  leading  to  ROC  curves  that  directly  compare 
performance  across  a  large  range  of  Pf  values. 

It  is  important  to  note  that  the  threshold  T  is  a  function  of  k.  This  implies  that  there  is 
a  different  threshold  for  each  of  the  hypotheses.  The  hypothesis  dependent  threshold  from 
equation  (4.15)  is  used  for  the  analysis  of  the  simulated  data.  In  the  next  section,  the  Bayes 
risk  is  derived  in  a  different  method,  creating  an  hypothesis  independent  threshold. 

4.2.5  Hypothesis  Independent  Threshold. 

Beginning  with  equation  (4.20),  the  goal  in  this  section  is  to  derive  a  detection 
algorithm  that  is  independent  of  k  in  the  threshold.  This  was  originally  proposed  to  help 
prove  that  the  prior  probability  of  the  null  hypothesis  does  not  affect  detection  performance. 
The  proof  of  the  null  hypothesis  prior  probability  is  covered  in  section  4. 2. 6. 2. 


(D(x,y)-B) 

x=l  y= 1 

Separating  equation  (4.16) 


Hk 

hk(x,y )  ^  In 

Ho 


cr 


Nx  Ny 

6  t  Z  h-(x,y) 

x=l  y=l 

2  cr 


(4.16) 


S 


,v=l  y=  I 


cr 

Nx  Ny 


hk(x,y ) 


Wt  = 


°",n(S) 


+ 


o  z  Z  h\x,y) 

x=ly=l 


(4.17) 


(4.18) 


6  2  cr 

Wu  is  effectively  a  weighting  term  that  modifies  the  SNR  value  based  on  the  likelihood 
of  the  hypothesis  occurring.  A  weighting  term  is  needed  for  each  hypothesis  k,  and  it  can 
be  precomputed  for  each  set  of  priors  and  PSF  hypotheses  used.  The  computation  of  these 
terms  does  not  add  significant  processing  time  compared  to  an  ECEP  test.  The  second 
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term  in  the  weighting  factor  is  what  is  traditionally  defined  as  the  threshold  term  T.  T  is 
generally  the  threshold  value  that  sets  the  false  alarm  rate  for  the  hypothesis  test.  This  is 
set  to  6  for  the  currently  implemented  SST  BHT,  but  it  is  modified  to  6.2212  to  keep  a 
constant  false  alarm  rate  in  the  ECEP  MHT  presented  in  [24].  More  discussion  regarding 
setting  the  desired  Pf  rates  is  covered  in  section  4.3.3. 

In  the  new  MHT  algorithm,  both  SNR  and  Wk  have  a  dependence  on  k,  and  are  moved 
to  the  left  side  of  the  equation.  Selecting  the  hypothesis  that  gives  the  largest  value  indicates 
which  hypothesis  is  most  likely  to  have  occurred.  This  hypothesis  is  then  compared  against 
the  null  hypothesis  to  make  a  final  detection  decision. 


Hk 

max  [SNR*  -  Wk]  ^  0  +  r  (4.19) 

k  Ho 

The  algorithm  is  theoretically  compared  against  a  threshold  of  zero,  but  a  threshold 
adjustment  term,  r,  is  added  to  achieve  different  Pf  values.  Selecting  a  specific  r  value 
allows  for  the  algorithm  to  operate  at  the  desired  false  alarm  rate.  Comparing  against  the 
proposed  ECEP  test  as  shown  below,  from  [24],  the  similarities  can  be  quickly  observed. 

Hk 

max  [SNR]  ^  T  (4.20) 

k  H0 

One  noticeable  benefit  is  that  the  equal  cost  portion  of  the  algorithm  still  preserves  the 
ability  to  select  a  hypothesis  k.  In  [67],  improved  detection  performance  was  observed,  but 
the  algorithm  could  only  make  a  binary  decision  between  a  space  object  being  present  or 
not.  The  sub-pixel  position  information  available  from  the  ECEP  and  ECUP  algorithms  is 
shown  to  be  useful  in  increasing  tracking  accuracy  by  Sligar  [69].  Another  difference  is 
the  additional  calculation  required  for  the  weighting  terms.  As  mentioned  previously,  W* 
does  not  depend  on  the  full  frame  data,  and  can  be  computed  once  per  frame  to  get  accurate 
background  standard  deviation  and  updated  priors  if  desired. 
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As  shown,  the  derived  ECUP  algorithm  is  able  to  be  reduced  to  a  sufficient  statistic, 
and  is  not  too  computationally  complex.  These  factors  are  all  positive  indicators  for  a 
successful  algorithm.  Next,  the  three  methods  for  determining  realistic  and  accurate  prior 
probabilities  are  investigated. 

4.2.6  Investigating  Null  Hypothesis  Probability. 

Two  methods  are  used  to  investigate  the  potential  values  for  the  null  hypothesis.  It 
is  theorized  that  the  value  selected  for  the  null  hypothesis  value  does  not  change  the 
capabilities  of  the  detection  algorithm.  If  this  is  shown  to  be  the  case,  there  is  no  need 
to  determine  an  accurate  estimate  for  kq. 

The  first  method  used  is  to  perform  the  detection  algorithm  on  simulated  data  with 
several  different  potential  nQ  values.  If  Pd  does  not  change  significantly  across  the  tested 
values,  it  can  be  assumed  that  the  value  for  no  has  no  impact.  Alternatively,  an  analytic 
approach  is  used  based  on  the  hypothesis  independent  MHT  proposed  in  equation  (4.19). 
These  approaches  are  covered  in  the  next  two  sections. 

4.2.6. 1  Simulation. 

After  assigning  alternate  hypotheses  prior  probabilities,  the  probability  for  null 
hypothesis  H0  needs  to  be  determined.  Accurately  determining  this  probability  empirically 
from  the  collected  data  is  difficult.  This  is  due  to  the  large  number  of  stars,  satellites,  debris, 
and  other  space  objects,  and  their  spatial  distribution.  Instead  of  attempting  to  estimate 
a  value  for  n0,  this  paper  shows  that  its  value  does  not  significantly  change  detection 
performance.  The  probability  of  an  RSO  being  present  in  any  pixel  being  tested  is  7Tq. 
Figure  4.3  shows  the  effect  that  implementing  the  detection  algorithm  from  equation  (4.20) 
with  no  values  ranging  from  0.05  to  0.95  has  on  computed  Pd.  Each  nQ  test  computed 
over  100  Monte  Carlo  iterations  to  determine  its  statistics.  Each  iteration  consists  of  100 
independent  RSOs  to  determine  the  Pd  value  at  the  specified  Pf  of  1  e  -  10. 
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Figure  4.3:  Analysis  of  P(i  for  varying  /r0  prior  probabilities  for  Pf  =  \e  -  10. 


Figure  4.3  shows  that  there  is  no  significant  change  in  Pd  for  a  wide  variety  of  no 
values.  In  addition,  investigating  the  means  and  standard  deviations  shows  that  all  n0  is 
well  within  one  standard  deviation.  As  a  result,  it  can  be  assumed  that  the  impact  of 
changing  the  null  hypothesis  prior  probability  no  is  negligible.  In  simulation,  the  value 
for  n o  is  chosen  to  be  0.10  to  match  the  prior  probability  in  the  ECEP  algorithm,  which 
assigned  all  priors  to  be  The  remaining  1  -  n0  probability  is  assigned  to  the  alternate 
hypotheses  with  the  proportions  described  in  Table  4.1. 

4.2.6.2  Analytic. 

All  the  methods  for  determining  prior  probabilities  for  the  null  hypotheses  described 
in  section  4.2.3  divide  100  percent  of  the  decision  space  between  the  M  hypothesis.  This 
forces  the  sum  of  n\  to  nM  to  be  one.  This  does  not  account  for  the  fact  that  the  total 
probability  of  all  the  hypotheses  must  sum  to  one,  including  n0.  To  accomplish  this,  the 
calculated  nk  are  adjusted  with  the  following  equation: 
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K  =  (1  -  n0 )nk. 


(4.21) 


If  there  is  no  analytical  or  empirical  method  of  selecting  a  value  for  no,  there  are  three 
potential  alternate  solutions.  One  would  be  to  use  the  same  assignment  as  ECEP,  nQ  =  0.10, 
since  it  is  the  method  being  compared  against.  Another  option  is  to  investigate  the  actual 
value  of  no,  but  as  mentioned  earlier,  this  may  not  be  possible.  Alternatively,  this  section 
demonstrates  analytically  that  the  effect  of  changing  the  null  hypothesis  prior  probability 
does  not  alter  detection  performance.  To  confirm  this  theory,  a  constant  C  is  introduced.  C 
acts  multiplicatively  with  n0,  to  adjust  the  null  hypothesis  value. 


7Tq  —  Cno 

n'k  =  (  \  -  Cn0)nk  (4.22) 

Proposing  different  values  for  no  can  be  accomplished  by  choosing  two  distinct 
constants,  C\  and  C2.  Computing  the  difference  in  weighting,  AW,,  gives  insight  into  the 
impact  of  changing  nQ. 


AW,  =  W,(C2)  -  W,(C,)  (4.23) 

W,(C)  is  the  weighting  term  from  equation  (4.18)  substituted  with  updated  priors  from 
equation  (4.22).  Combining  and  reducing  with  constants  C\  and  C2  gives  an  expression  for 
the  change  in  the  kth  weighting  term  due  to  a  change  in  the  null  hypothesis  prior  probability. 


AW,  = 


—In  [  C*“  ) 

6  \(1  -  C2nQ)nk) 


—Inf  C'm  ) 

6  \(1 -Ci7r0)7r,/ 


(4.24) 


AW,  =  -In 
6 


C2(l  -  Ci7Tq)\ 
Cxi  1  -  C2no)l 


(4.25) 
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Looking  at  equation  (4.25),  AW/.  is  only  a  function  of  the  constants  and  the  original 
ttq.  This  signifies  that  A Wk  is  independent  of  k,  and  it  can  be  moved  and  grouped  into 
the  threshold  r.  The  constant  adjustment  to  the  threshold  will  change  the  probability  of 
false  alarm.  The  false  alarm  probability  can  then  be  adjusted  to  match  the  Pf  rate  of  the 
algorithm  being  compared  against.  This  analysis  demonstrates  that  the  assignment  for  the 
value  of  7r0  does  not  change  the  detection  performance  of  the  algorithm. 

Combining  the  simulation  study  on  no  sensitivity  along  with  the  analytic  derivation 
above  provides  good  evidence  that  the  choice  of  null  hypothesis  prior  probability  does  not 
impact  the  detection  capabilities  of  the  algorithm. 

4.3  Experiment  Descriptions 

In  this  section,  the  simulations  and  experiments  used  to  test  the  algorithms  developed 
in  section  4.2.4  and  4.2.5  are  described. 

4.3.1  Simulation  Description. 

To  test  the  newly  developed  detection  algorithm,  simulated  RSOs  are  analyzed.  The 
primary  reason  simulated  data  is  used  is  because  a  large  number  of  space  objects  uniformly 
distributed  across  the  pixels  being  tested  are  needed  to  test  the  algorithm.  It  is  difficult  to 
assess  algorithm  performance  by  analyzing  collected  telescope  data  of  a  single  RSO.  It 
is  possible  that  the  object  is  located  in  a  specific  sub-pixel  position  on  the  CCD  and  does 
not  move  throughout  the  collection  period.  In  this  case,  the  object  may  not  represent  how 
other  randomly  located  PSFs  appear  in  the  data.  Alternatively,  another  potential  method 
is  to  analyze  a  large  image  of  the  sky,  including  stars,  to  gather  a  variety  of  locations  and 
intensity  levels.  A  drawback  to  this  method  is  that  many  of  the  bright  objects  will  be  easily 
detectable  to  both  algorithms,  therefore  it  is  necessary  to  define  a  range  of  intensities  where 
the  algorithms  perform  differently. 

In  this  simulation,  it  is  assumed  that  the  RSOs  occur  with  the  prior  probability  that  is 
present  in  an  actual  data  collection.  The  ECUP  algorithm  has  the  maximum  performance 
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advantage  compared  to  ECEP  if  the  true  priors  exactly  match  the  priors  determined  through 
the  correlation  method.  The  priors  may  differ  in  an  actual  data  collection  if  the  assumption 
of  a  uniform  distribution  of  objects  does  not  hold,  or  other  factors  not  included  in  the 
analysis  affect  the  probability.  In  this  analysis,  the  ideal  performance,  which  does  not 
include  those  factors,  is  tested.  To  give  the  simulated  data  the  desired  prior  probabilities, 
a  function  is  created  to  randomly  select  one  of  the  hypotheses  from  the  nine  alternate 
hypotheses  with  the  assigned  weighting  probabilities  The  optical  model  described  in 
section  4.2.1,  h(x,y),  is  used  to  create  the  9  candidate  PSFs.  As  mentioned,  this  model 
includes  several  optical  effects  including  a  long  exposure  atmosphere  model  [34],  optical 
aberrations,  and  spatial  aliasing  of  the  CCD  pixels.  The  data  model  in  equation  (4.8)  is 
used  to  combine  the  PSF  with  the  background  and  noise  model. 

Two  types  of  simulated  data  are  needed,  one  with  a  RSO  present  and  one  without 
a  RSO.  Both  data  sets  are  needed  to  create  the  ROC  curves  in  section  4.4.2.  For  this 
simulation,  an  object  photocount  of  Q  -  300  and  background  level  of  B  -  400  are 
used.  In  the  case  where  no  object  is  present,  B  is  still  200,  but  6  is  set  to  zero.  To  add 
noise,  the  PSF  and  background  is  set  to  the  mean  of  a  Poisson  random  variable.  Using  a 
Poisson  noise  assumption  implies  that  the  signal  is  dominated  by  photon  counting  noise. 
Other  factors,  including  dark  current  and  read  noise,  are  not  significant  contributors.  The 
detection  algorithm  assumes  that  the  received  noise  is  Gaussian,  as  described  in  section 
4.2.4.  Figure  4.4  shows  a  flow  diagram  of  the  steps  used  to  generate  simulated  data  with 
and  without  an  RSO  present. 

4.3.2  Collected  telescope  data. 

The  second  experimental  method  used  is  collected  SST  data.  This  experiment  uses 
the  same  data  set  described  in  Chapter  3,  but  focuses  on  a  different  portion  of  the  data 
present  in  the  images.  This  data  was  originally  collected  as  an  experiment  to  test  how  faint 
of  an  object  can  be  detected  by  different  detection  algorithms.  The  SST  is  programmed  to 
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Figure  4.4:  Flow  diagram  of  process  to  (A)  simulate  data  with  an  RSO  present  and  (B) 
simulate  data  with  no  RSO  present.  Both  data  sets  are  needed  to  generate  ROC  curves. 


track  a  communication  satellite,  ANIK-F1,  as  it  enters  eclipse.  The  resulting  intensity  loss 
causes  the  object  to  go  from  easily  detectable  to  difficult  to  detect.  This  is  used  to  test  the 
performance  of  the  UCEP  algorithm. 

In  addition  to  the  satellite  of  interest  in  each  image,  there  are  thousands  of  other 
objects,  including  stars  and  potentially  RSOs.  These  stars  and  other  space  objects  give 
a  large  quantity  of  varying  intensity  objects  to  test  algorithms  against.  By  processing  a 
large  portion  of  an  image  and  totaling  the  number  of  detections,  a  metric  of  performance 
can  be  determined  for  each  algorithm.  Comparing  the  total  number  of  detections  with  each 
algorithm  for  the  same  false  alarm  rate  gives  one  method  to  compare  the  algorithms.  All 
types  of  objects  are  treated  similarly  in  this  experiment,  since  both  stars  and  RSOs  appear 
as  unresolved  point  source  objects  to  the  SST. 

Each  frame  of  collected  SST  data  consists  of  6144x4096  pixels,  where  each  pixel 
contains  2x2  binned  1 5//m  square  pixels.  The  data  analyzed  in  this  paper  was  collected  on 
three  nights  13-15  March  2012,  referred  to  as  night  073-075  respectively  in  the  following 
discussions.  As  mentioned  in  the  optical  model  description,  each  image  is  a  long  exposure 
collection.  This  implies  that  a  long  exposure  atmosphere  is  present  in  the  data. 
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Within  each  image  there  are  areas  where  CCD  arrays  are  aligned,  and  other  areas 
of  imperfection.  To  remove  these  edge  effects,  a  subset  of  each  frame  is  analyzed.  An 
area  consisting  of  1024x1024  pixels  is  selected  from  each  image.  The  same  pixels  are 
consistently  used  from  each  frame  as  the  data  is  processed  by  the  algorithm.  Testing  on 
each  pixel  is  done  by  including  a  16x16  window  around  the  pixel  being  tested.  Including 
the  window  gives  enough  samples  to  calculate  accurate  background  statistics  as  well  as 
enough  pixels  to  capture  both  the  PSF  model  and  any  potential  space  objects.  In  the  next 
section,  the  false  alarm  probabilities  are  investigated. 

4.3.3  Setting  False  Alarm  Probability  for  Collected  Data. 

It  is  important  to  ensure  that  both  the  ECEP  and  ECUP  algorithms  have  the  same 
probability  of  false  alarm,  Pf.  The  weighting  term,  W*,  is  effectively  changing  the  threshold 
and  therefore  the  false  alarm  probability.  Since  the  weighting  term  can  be  different  for  each 
hypothesis,  the  false  alarm  probability  for  the  entire  test  needs  to  account  for  all  potential 
hypotheses. 

To  accurately  calculate  Pf,  two  assumptions  are  made.  The  first  is  that  SNR  statistics 
calculated  with  equation  (4.17)  are  Gaussian,  which  follows  logically  from  the  assumption 
that  the  received  noise  is  Gaussian.  The  second  assumption  is  that  the  probabilities  of 
false  alarm  for  each  hypothesis  are  independent,  and  can  be  calculated  separately.  This 
assumption  is  applied  to  both  the  ECEP  and  ECUP  algorithms.  This  assumption  may 
overestimate  the  expected  Pf  if  there  is  any  dependence  between  hypotheses,  but  will  do 
so  for  both  the  ECEP  and  ECUP  algorithms.  For  the  ECEP  test,  a  base  threshold  of  y  =  6 
is  used  to  match  previous  algorithms  and  the  SST  program’s  threshold.  The  threshold  for 
ECUP  is  then  adjusted  to  match  the  Pf  produced  by  the  threshold  of  y  =  6. 

To  mathematically  describe  the  method  of  calculating  Pf,  the  term  N(pi,  cr )  is  defined 
to  represent  a  Gaussian  distribution  with  mean  /i  and  standard  deviation  cr.  As  described 
previously,  the  SNR  in  the  ECEP  algorithm  is  normalized  to  be  a  zero  mean  unit  variance 
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Gaussian  random  variable.  Since  the  ECUP  algorithm  is  not  normalized  by  <rk,  SNRA  is  no 
longer  zero  mean  unit  variance  Gaussian.  SNRA  has  a  mean  of  Wk  and  standard  deviation 
of  crk  for  each  hypothesis. 


ECEP:  X  ~  JV( 0, 1) 

ECUP:  Xk~N(Wk,o-k ) 

Assuming  independence  between  the  hypothesis  tests,  the  total  false  alarm  can  be 
calculated  by  summing  the  false  alarm  probability  from  each  alternate  hypothesis. 

M 

Pf  =  J]l-p(X<T)  (4.26) 

k=  1 

p(X  <  t)  is  the  Cumulative  Distribution  Function  (CDF),  and  X  is  a  random  instance 
drawn  from  the  Gaussian  distribution  described  by  either  the  ECUP  or  ECEP  test.  For  the 
ECEP  test,  a  value  of  6  is  used  for  r.  To  ensure  the  tests  have  the  same  probability  of 
false  alarm,  a  r  needs  to  be  found  that  gives  equal  Pj  between  the  tests.  Next,  the  results 
for  the  calculated  prior  probabilities  and  the  number  of  binary  detections  are  presented  and 
analyzed. 

4.4  Results 

This  section  will  investigate  the  results  from  both  the  simulated  experiment  and  the 
data  collected  from  the  SST. 

4.4.1  Decision  Space  Results. 

The  first  important  result  is  the  calculated  values  for  the  prior  probabilities  nk.  To 
calculate  these  probabilities,  the  total  number  of  each  hypothesis  present  in  the  hypothesis 
assignment  matrix  PL  is  divided  by  the  total  number  of  sub-pixel  positions  tested.  It  is 
helpful  to  define  the  mathematical  notation  of  an  Iverson  Bracket  [70],  where  [•]  gives 
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a  value  of  1  is  true  or  zero  if  false  the  following  equation  gives  the  value  for  the  prior 
probabilities: 


nk 


Na  Na 

XX 


<2=1  OJ=  1 


ma,io)  =  Hk\ 

NaNw 


(4.27) 


Na  and  Nw  are  the  size  of  the  window  being  investigated  in  a  and  to.  Equation  (4.27) 


is  applied  to  each  hypothesis  assignment  matrix.  The  three  methods  calculating  the  prior 


probabilities  are  shown  in  Figure  4.5. 
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Figure  4.5:  Comparison  of  three  proposed  methods:  distance,  correlation,  and  empirical 
based  methods  of  assigning  sub-pixel  positions. 
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Figure  4.5  shows  the  prior  probabilities  for  the  three  methods  of  segmenting  the 
decision  space.  The  first  detail  to  note  is  that  all  three  methods  give  similar  prior  probability 
values,  with  a  maximum  difference  of  approximately  10  percent.  The  general  agreement 
between  the  priors  gives  a  good  indication  that  the  methods  accurately  represent  the  true 
probabilities.  The  hypothesis  at  the  center  of  the  pixel,  H5,  has  the  largest  probability  of 
occurring  in  all  methods.  The  prior  value  ranges  between  25  to  30  percent.  Using  a  BHT, 
only  one  hypothesis  is  considered,  and  it  is  that  the  object  is  located  directly  in  the  center  of 
a  pixel.  This  analysis  shows  that  using  a  BHT  would  only  match  well  with  25  to  30  percent 
of  potential  objects,  decreasing  the  ability  to  detect  space  objects. 

The  methods  vary  slightly  in  how  they  assign  probability  between  the  sides  and 
comers.  Distance  and  correlation  metrics  give  higher  probability  to  the  sides,  while  the 
empirical  method  favors  the  corners.  The  empirical  priors  for  the  sides  and  comers  are 
closer,  while  the  other  methods  have  a  larger  separation  between  them.  These  differences 
can  be  due  to  things  not  accounted  for  in  the  model  that  occur  in  measured  data,  such  as 
non-uniform  response  across  the  CCD  detector  areas  within  pixel.  In  the  empirical  method, 
there  is  noise  present  in  the  data  which  is  not  present  with  the  other  methods. 

Considering  the  three  methods  and  their  resulting  priors,  the  empirical  method  of 
determining  tt&  values  is  used  for  analyzing  collected  data.  This  is  because  it  is  the 
closest  match  to  how  the  ECUP  detection  algorithm  actually  operates.  Additionally,  it 
includes  noise  and  other  details  that  are  not  present  in  the  other  methods.  The  distance  and 
correlation  metrics  also  provide  confidence  that  the  empirical  method  is  an  accurate  way  to 
define  prior  probability  since  they  generally  agree.  Table  4. 1  summarizes  important  details 
about  the  physical  locations  of  the  hypotheses  within  a  30x30  pixel,  prior  probabilities  with 
the  empirical  metric,  and  the  resulting  weighting  values. 

As  Table  4.1  shows,  probability  is  lower  in  the  comers  and  sides,  with  the  additional 
probability  shifted  to  H5,  or  the  center  pixel  hypothesis.  This  implies  that  it  is  more  likely 
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Table  4.1:  The  sub-pixel  locations  of  the  M  hypotheses  (join )  along  with  the  calculated 
priors  nk  and  weighting  values  Wk  based  on  the  empirical  method. 


Hi 

h2 

h3 

h4 

h5 

h6 

h7 

h8 

h9 

a,  co 

-15,-15 

-15,0 

-15,15 

0,-15 

0,0 

0,15 

15,-15 

15,0 

15,15 

TCk 

0.10 

0.08 

0.10 

0.09 

0.26 

0.09 

0.10 

0.08 

0.10 

Wk 

1.78 

2.57 

1.78 

2.56 

3.75 

2.56 

1.78 

2.56 

1.78 

for  captured  data  to  closely  match  the  PSF  generated,  assuming  the  object  is  in  the  center 
of  the  pixel.  Based  on  the  empirical  method,  H5  hypothesis  has  approximately  a  26  percent 
chance  of  occurring.  The  weighting  values,  Wk,  from  equation  (4.18)  which  are  calculated 
based  on  nk,  also  provide  interesting  results.  All  of  the  weighting  values  are  positive,  and 
Wk  is  subtracted  from  the  computed  SNR.  The  newly  computed  SNR  compared  against  the 
threshold  will  be  lower  than  the  ECEP  test.  This  is  counteracted  by  the  new  threshold  r 
computed  to  ensure  equal  Pf  rates.  Although  the  weighting  term  is  a  linear  effect  on  the 
SNR,  calculating  r  is  not  linear.  Considering  these  two  competing  effects,  it  is  difficult 
to  analyze  how  this  will  increase  or  decrease  the  number  of  RSOs  detected.  This  is  why 
collected  telescope  data  detection  performance  is  analyzed. 

4.4.2  Processing  Simulated  Data  with  MHT  Algorithms. 

In  this  section,  the  detection  algorithm  is  applied  to  simulated  RSO  data.  As  has 
been  done  in  previous  chapters,  using  a  ROC  curve  gives  insight  to  the  detection  behavior 
over  a  large  range  of  potential  Pf  values.  For  the  simulated  data,  the  prior  probabilities 
determined  through  the  correlation  metric  were  used.  This  is  because  the  empirical  method 
is  not  realistic,  because  there  is  no  collected  data  to  process.  The  following  equation  is  used 
to  generate  ROC  curves: 
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(4.28) 


max  [SNRa.  -  Wk\  %  Y  +  r. 

k  H0 

t  is  varied  to  give  Pc\  and  Pf  pairs.  The  ability  to  correctly  identify  the  sub-pixel 
hypothesis  is  not  reported.  If  data  containing  a  RSO  is  processed  and  any  of  the  alternate 
hypotheses  is  selected,  it  is  counted  as  a  correct  detection.  Figure  4.6  shows  an  instance 
of  ROC  curves  for  both  the  ECEP  and  ECUP  algorithms.  To  generate  these  ROC  curves, 
100  independent  simulated  images  are  analyzed.  For  every  image,  the  SNR  value  and 
weighting  factor  wk  is  calculated  for  all  hypotheses.  The  maximum  SNR  value  is  selected 
and  compared  against  the  threshold  T. 


Figure  4.6:  A  single  instance  of  a  ROC  curve  for  both  the  ECEP  and  ECUP  algorithms. 
Simulated  RSO  and  data  with  no  objects  present  are  needed  to  generate  probability  of 
detection  and  false  alarm. 
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The  ROC  curves  in  figure  4.6  show  up  to  approximately  10  percent  of  ECUP 
over  ECEP.  The  algorithms  perform  similarly  for  relatively  high  and  low  false  alarm 
probabilities.  A  typical  detection  algorithm  would  not  operate  in  these  regimes.  The  two 
curves  in  figure  4.6  represent  a  single  instance  of  a  simulation.  To  ensure  that  the  ECUP 
algorithm  detects  more  objects  than  the  ECEP  generally,  the  statistics  of  the  ROC  curve  are 
investigated. 

Defining  A pd(Pf)  as  the  difference  between  Pd  for  the  ECUP  and  ECEP  algorithms 
gives  the  detection  difference  for  the  same  false  alarm  probability. 

Ap/Pf)  =  ECUP  -  ECEP  (4.29) 

Figure  4.7  shows  a  statistical  comparison  of  detection  performance  between  ECEP 
and  ECUP  algorithms.  The  Pd  values  in  this  figure  were  calculated  with  Pj  =  \e  -  10. 

As  figure  4.7  shows,  there  is  a  small  but  statistically  significant  difference  in  detection, 
A Pd,  between  the  ECEP  and  ECUP  algorithms.  The  red  dashed  line  indicates  the  mean 
difference  across  100  Monte  Carlo  simulations.  The  improvement  in  this  simulation 
is  approximately  5  percent.  There  are  instances  where  APd  is  negative,  or  the  ECEP 
outperforms  the  ECUP  algorithm,  but  on  average  this  is  not  the  case.  A  5  percent 
improvement  across  a  large  number  of  RSO  could  still  be  a  significant  number  of  objects. 
In  the  next  section,  the  algorithm  is  used  to  process  collected  SST  data. 

4.4.3  Processing  SST  Data  with  MHT  Algorithms. 

After  assigning  priors,  the  ECUP  algorithm  can  be  implemented  on  collected  data. 
To  accomplish  this,  the  data  described  in  section  4.3.2  is  analyzed.  This  data  contains 
many  objects  of  different  intensity  levels,  including  many  much  greater  than  the  current 
threshold.  These  objects  will  be  easily  detectable  with  both  the  ECEP  and  ECUP  algorithm. 
The  difference  in  detection  is  the  RSOs  that  are  close  to  the  detection  threshold.  Objects 
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Figure  4.7:  Statistical  comparison  of  ECEP  and  ECETP  algorithms  over  100  monte  carlo 
iterations.  Pd  values  for  both  algorithms  are  found  for  a  false  alarm  probability  of 
Pf  =  le  -  10. 


much  brighter  than  these  will  easily  be  detected  with  both  algorithms.  Space  objects  much 
dimmer  than  the  threshold  will  not  be  detected  by  either  algorithm. 

There  are  several  potential  methods  for  reporting  the  change  in  detection  performance. 
Using  a  probability  of  detection,  Pd,  at  a  specified  false  alarm  rate  Pf  is  a  commonly  used 
metric.  Another  method  is  reporting  Pd  across  a  large  range  of  potential  Pf  values,  also 
known  as  a  ROC  curve,  typically  provides  insight.  In  this  chapter,  the  number  of  additional 
binary  detections  at  a  specified  false  alarm  is  reported.  The  difference  in  the  number  of 
detected  objects  between  ECEP  and  ECUP  is  defined  as  A0: 
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A0  =  ECUP  -  ECEP. 


(4.30) 


The  total  number  of  binary  detections  includes  all  alternate  hypotheses  where  an  object 
is  considered  present,  H\  through  H9.  Using  the  calculated  Pf  value  of  8.5e-09,  and  r  =  6 
for  ECEP  and  r  =  0.98  for  ECUP,  the  SST  data  is  analyzed.  Figure  4.8  demonstrates  a 
small  subset  of  a  frame  of  data  processed  with  both  the  ECEP  and  ECUP  algorithms.  In 
the  example  shown,  one  additional  binary  detection  is  observed. 


Figure  4.8:  A  sample  subset  of  SST  data  processed  with  (A)  the  ECEP  algorithm  and  (B) 
the  ECUP  algorithm.  In  this  example,  one  additional  detection  is  observed. 


Figure  4.9  demonstrates  the  difference  in  the  number  of  binary  detections  between 
ECUP  and  ECEP. 

A0  ranges  from  61  additional  binary  detection  to  4  less  detections  on  the  days  and 
frames  processed.  Across  all  three  nights  processed  there  were,  on  average,  more  detections 
with  ECUP.  For  each  frame,  the  total  number  of  detections  changes  as  objects  enter  and  exit 
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Figure  4.9:  Difference  in  number  of  binary  detections,  A0,  over  32  sample  frames  from 
data  collected  on  nights  (A)  073,  (B)  074,  and  (C)  075. 


the  frame,  noise  spikes  occur,  and  space  objects  change  intensity.  A0  stays  fairly  consistent 
throughout  these  changes. 

There  is  an  obvious  increase  in  binary  detections  present  in  the  sample  data  analyzed, 
but  it  is  important  to  ensure  that  this  instance  is  not  a  random  occurrence.  One  method 
is  to  perform  a  significance  analysis  on  the  mean.  Looking  at  a  paired  T  statistic  test  can 
indicate  if  there  is  truly  a  difference  in  means,  and  which  is  greater  [71].  Night  073  has  a 
paired  T  test  statistic  value  of  9.07.  This  signifies  it  is  highly  likely  that  ECUP  is  detecting 
more  objects  on  average  than  ECEP.  Large  T  statistic  values  also  occur  on  nights  074  and 
075,  15.62  and  9.49  respectively. 

Additional  information  can  be  gleaned  from  investigating  the  distribution  of  A0  data. 
Figure  4.10  shows  a  histogram  for  A0  across  all  three  nights  of  the  SST  data  processed. 

The  distribution  of  the  additional  number  of  detections  appears  Gaussian.  Across  all 
of  the  processed  data,  Aa  has  a  mean  of  22.73  and  a  standard  deviation  of  12.65.  Table  4.2 
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Figure  4.10:  Combined  histogram  for  A„,  the  number  of  additional  binary  detections  with 
ECUP  algorithm.  All  three  nights  analyzed  073-075  are  included. 


shows  the  average  number  of  additional  detections  and  the  standard  deviation  of  the  ECUP 
algorithm  for  individual  nights. 


Table  4.2:  Average  number  of  additional  binary  detections  by  ECUP,  A0,  and  standard 
deviation,  crAo,  for  three  nights  13-15  March  2012. 


Night  073  Night  074  Night  075 

A0  26.13  22.97  19.10 

mAo  16.30  8.32  11.38 
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One  noticeable  result  from  Table  4.2  is  that  the  largest  A0  occurred  on  night  073. 
This  was  the  night  used  to  generate  the  priors  with  the  empirical  method.  This  signifies 
two  interesting  details.  The  first  is  that  the  empirical  priors  do  not  necessarily  need  to  be 
recomputed  frequently.  The  empirical  priors  are  able  to  achieve  an  improvement  for  three 
consecutive  nights  without  adjustment.  The  values  do  decrease  slightly  between  nights. 
Several  factors  including  atmospheric  conditions  may  cause  the  priors  to  change  slightly 
between  collections.  There  may  be  an  optimal  update  rate  for  priors  that  takes  these  factors 
into  account,  but  this  is  not  further  investigated  in  this  paper. 

The  subset  of  frames  analyzed  for  this  chapter  are  not  the  full  image  frames.  Increasing 
the  number  of  pixels  being  tested  would  also  increase  A0.  A  total  of  24  of  the  1024x1024 
frames  are  in  the  full  frame  data.  Applying  the  averages  determined  in  the  small  frame 
would  give  approximately  458-627  additional  binary  detections  on  average. 

4.5  Full  Frame  Implementation 

An  important  consideration  for  implementing  an  algorithm  is  how  long  it  takes  to 
process  each  frame  of  data.  If  the  algorithm  can  not  be  completed  in  real  time,  or  the  time 
between  when  each  image  is  captured,  the  output  of  the  algorithm  would  be  constantly 
falling  behind  the  captured  data.  Some  catch  up  can  be  accomplished  when  data  collections 
are  not  being  performed.  This  limits  the  ability  to  adjust  collections  based  on  observations 
during  collection. 

In  section  3.5,  a  full  frame  implementation  of  the  UCEP  algorithm  is  discussed. 
Similar  methods  can  be  used  with  the  UCEP  algorithm.  The  main  difference  between 
full  frame  implementation  of  UCEP  and  ECUP  algorithms  is  exponential  function.  Unlike 
the  UCEP  algorithm,  which  is  not  a  linear  algorithm,  the  ECUP  algorithm  is  able  to  be 
reduced  to  SNR  calculation,  with  an  additional  weighting  term.  This  weighting  term,  W*, 
can  be  precomputed  based  on  the  background  statistics,  PSF  being  considered  fk(x,y),  and 
the  assigned  priors.  These  terms  do  not  add  a  significant  processing  requirement. 
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4.6  Conclusions 


In  this  chapter,  the  assignments  of  a  priori  probabilities  in  a  MHT  are  investigated. 
This  is  achieved  through  developing  an  algorithm  and  testing  it  against  simulated  and 
collected  SST  data. 

To  analyze  the  simulated  data,  a  correlation  metric  between  the  simulated  RSO  data 
D(x,y )  and  the  modeled  sub-pixel  positions  Ta^(x,y),  is  used  to  assign  prior  probabilities 
for  the  M  hypotheses.  It  is  shown  with  a  sensitivity  analysis  that  the  null  hypothesis 
does  not  significantly  impact  the  detection  ability,  and  a  value  of  0.10  is  used  to  match 
previously  developed  algorithms.  These  prior  probabilities  are  implemented  in  an  ECUP 
MHT.  The  derived  algorithm  is  used  to  analyze  simulated  RSOs.  The  algorithm  used  for 
this  research  contains  a  hypothesis  dependent  threshold.  ROC  curve  analysis  shows  that  the 
ECUP  algorithm  has  a  greater  Pd  than  the  ECEP  algorithm  for  a  large  range  of  Pf  values. 
At  a  relevant  false  alarm  level  of  Pf  =  10-10,  ECUP  shows  an  average  improvement  of 
approximately  5  percent  when  tested  over  100  monte  carlo  simulations. 

When  tested  against  collected  SST  data,  this  research  demonstrates  that  for  certain 
applications  it  may  be  beneficial  to  use  an  ECUP  MHT  algorithm.  It  is  shown  that 
with  high  statistical  confidence  the  ECUP  algorithm  on  average  has  more  detections  than 
an  ECEP  algorithm  performed  on  the  SST  data.  For  the  small  subset  of  the  frames 
analyzed  approximately  19-26  additional  detections,  which  translates  to  458-627  additional 
detections  on  average  for  full  frame  data.  These  detections  are  most  likely  threshold 
objects  that  are  difficult  for  currently  proposed  detection  algorithms  to  detect.  Finding  these 
threshold  objects  will  provide  additional  information  for  tracking  and  characterization,  in 
turn  increasing  SDA  capabilities  with  the  optical  system. 

In  addition,  a  method  of  using  empirical  observations  of  representative  data  is 
developed  to  form  a  model  for  a  priori  probabilities.  Distance  and  correlation  metrics  were 
also  considered  and  result  in  similar  distribution  of  prior  probabilities. 
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The  ECUP  algorithm  assumes  Gaussian  received  noise  in  this  paper.  In  the  future, 
removing  the  Gaussian  assumption  may  detect  additional  space  objects,  if  the  received 
noise  is  not  truly  Gaussian.  Additional  methods  may  also  be  combined  to  potentially 
improve  detection  performance  even  further.  For  example,  using  a  unequal-cost  and  an 
unequal-prior  probability  approach  may  improve  over  a  unequal-cost  equal-prior  test. 

In  the  ECEP  MHTs,  the  threshold  is  considered  constant  for  each  hypothesis.  As  a 
result,  the  modeled  intensity  is  effectively  different  for  each  hypothesis,  and  is  dependent 
on  the  distribution  of  the  PSF.  This  assumption  is  changed  in  this  paper  for  the  ECUP 
algorithm.  This  assumption  can  be  removed  in  a  future  ECEP  test,  and  may  provide 
interesting  research  opportunities  and  advancements. 

Another  consideration  is  how  including  unequal-priors  may  effect  the  ability  to 
translate  the  detected  objects  into  accurate  tracking  and  orbit  determination.  The  ECUP 
algorithm  takes  into  account  the  probability  of  the  alternate  hypotheses,  each  of  which 
corresponds  to  a  sub-pixel  position.  The  sub-pixel  information  gathered  from  performing 
a  ECUP  test  may  lead  to  improved  sub-pixel  position,  and  in  turn  more  accurate  tracking 
of  detected  RSOs. 
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V.  Characterizing  Point  Spread  Function  fluctuations  to  improve  Resident  Space 

Object  detection 


In  the  previous  chapters,  this  dissertation  examines  Bayes  Risk  algorithms  based  on 
the  assumption  that  the  received  data  contains  Gaussian  noise.  This  chapter  focuses  on 
developing  an  accurate  combined  noise  statistics  model  for  noise  present  in  a  collected 
telescope  image.  The  specific  goal  of  this  chapter  is  to  determine  a  more  accurate  model 
for  the  observed  noise  in  a  PSF.  By  using  this  new  model  and  the  corresponding  detection 
algorithms,  there  will  be  improved  detection  performance  in  the  optical  systems  being 
considered.  As  mentioned  previously,  the  PSF  is  of  interest  due  to  the  fact  that  all  objects 
of  interest  can  be  considered  point  source  objects. 

Building  a  detection  algorithm  based  on  a  more  accurate  or  complete  model  of  noise 
statistics  in  a  PSF  should  theoretically  result  in  an  increase  in  the  ability  to  detect  space 
objects  with  that  optical  system.  The  number  of  additionally  detected  objects  and  the 
complexity  of  the  resulting  algorithm  are  important  factors  in  determining  the  feasibility  of 
any  newly  proposed  detection  algorithm. 

This  chapter  is  organized  into  four  main  sections.  In  section  5.1,  relevant  background 
and  methodology  in  optical  modeling,  detection  theory,  atmospheric  effects,  and  relevant 
Probability  Density  Functions  (PDFs)  are  detailed.  Section  5.2  describes  the  simulation 
experiment  used  to  analyze  the  distributions  and  detection  performance.  Section  5.3 
discusses  the  results  and  analysis,  and  section  5.4  presents  the  conclusions  on  the  results 
reached. 

5.1  Methodology 

All  objects  of  interest  in  this  research,  due  to  their  relative  size  and  distance  from  the 
observation  plane,  appear  as  point  sources  to  the  telescope.  The  PSF  of  an  optical  system 
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is  the  response  of  the  system  to  an  impulse,  or  point  source.  In  the  case  of  RSO  detection, 
each  object  is  an  observed  PSF. 

The  noise  present  in  an  image  captured  through  a  telescope  is  a  combination  of  several 
potential  noise  sources.  Traditionally,  SDA  detection  algorithms  have  assumed  a  Central 
Limit  Theorem  justification  and  used  a  Gaussian  PDF  to  describe  the  distribution  of  the 
received  data.  This  implies  that  the  combination  of  all  of  the  noise  sources  present  in  the 
system  result  in  an  overall  Gaussian  PDF  within  the  received  data.  Multiple  research  efforts 
have  looked  at  creating  detection  algorithms  based  on  the  noise  in  the  received  data  not 
being  Gaussian  [22,  35].  In  these  cases  it  is  assumed  that  the  received  noise  is  dominated 
by  Poisson  noise,  and  other  sources  do  not  contribute  significantly. 

One  example  of  building  a  detector  based  on  Poisson  noise  is  when  the  dominant 
source  of  noise  in  the  received  data  is  assumed  to  be  the  random  arrival  time  of  photons 
[22].  The  author  states  this  can  be  achieved  through  a  very  low  noise  CCD  in  a  low 
intensity  collection  environment.  It  is  found  that  this  assumption  reduces  similarly  to 
Gaussian  if  there  is  a  flat  background.  In  other  words,  if  the  detection  algorithm  does 
not  have  a  known  star  map  for  use  by  the  detection  algorithm,  the  Poisson  and  Gaussian 
approaches  perform  similarly.  Similar  results,  also  assuming  a  Poisson  distribution,  have 
been  demonstrated  to  perform  comparably  in  collected  telescope  data  at  the  Air  Force 
Institute  of  Technology  (AFIT)  [35]. 

Unlike  to  these  two  non-Gaussian  methods,  this  chapter  will  account  for  two  of  the 
main  phenomena  present  in  data  collection,  and  combine  them  into  a  single  noise  model. 
These  phenomena  are  atmospheric  effects  and  photon  counting  noise.  It  is  well  understood 
that  noise  due  to  random  photon  arrival  times  follows  a  Poisson  distribution  [17].  The 
impact  of  the  atmosphere  and  the  fluctuations  it  causes  in  received  intensity  are  not  as 
clearly  defined. 
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5.1.1  Optical  Model. 

As  in  Chapter  3,  an  optical  model  of  the  telescope  is  needed  to  generate  the  PSFs 
used  in  the  detection  algorithm.  There  are  two  main  differences  in  the  optical  model 
presented  here:  one  is  the  method  for  simulating  the  atmospheric  effects  and  the  second 
is  the  exclusion  of  lens  aberrations.  This  implies  that  the  optical  system  being  modeled 
is  either  diffraction  limited  or  that  the  lens  aberrations  do  not  change  significantly  over 
the  time  scale  of  interest.  The  PSF  of  an  optical  system,  hopt(x,y),  can  be  found  by  the 
following  relation  [34]: 


hopt(x,y)  =  \T{P(m,n)}\2  (5.1) 

where  x,y  are  spatial  distance  pixel  coordinates  in  the  detector  plane,  and  m,n  are  pixel 
coordinates  in  the  pupil  plane.  P(m,  n )  is  a  pupil  function  that  mathematically  describes  the 
effect  of  the  pupil  on  incoming  light  and  T  is  a  two-dimensional  Fourier  transform.  In  a 
diffraction  limited  optical  system  consisting  of  a  perfect  lens  or  mirror,  the  pupil  function 
consists  of  only  the  geometry  of  the  pupil,  P(m,  n )  =  A  (in,  n),  and  contains  no  other  phase 
distortions.  A(m,  n)  is  an  amplitude  function  that  is  one  or  zero,  depending  on  if  light  is 
able  to  pass  through  the  pupil  at  a  specific  m,  n  pixel  location. 

In  a  non-ideal  imaging  system,  imperfections  in  the  lenses  or  mirrors  cause  phase 
distortions  to  any  light  passing  through  the  optics.  These  distortions  are  modeled  as  phase 
fluctuations  to  the  amplitude  function  [34]. 

P(m,  n)  =  A(m,  n)  exp  [j0o(m,  n)]  (5.2) 

For  this  simulation,  it  is  assumed  that  the  imaging  system  is  ideal.  This  assumption 
implies  that  the  phase  function  is  a  constant  function  of  spatial  coordinates  m  and  n.  This 
is  justifiable  due  to  the  fact  that  the  phase  distortions  caused  by  the  optics  are  generally 
constant  in  time  and  these  distortions  will  not  impact  the  peak  intensity  as  a  function  of 
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time.  The  lens  aberrations  change  the  shape,  distribution,  and  intensity  of  the  PSF,  but  do 
not  generally  change  significantly  over  the  observation  time. 

At  this  point,  the  optical  model  represents  a  field  entering  the  pupil  of  an  optical 
system,  but  has  not  accounted  for  the  effects  of  propagation  through  the  atmosphere.  To 
randomly  generate  a  phase  screen  with  statistics  that  accurately  represent  the  atmosphere, 
a  Zernike  method  is  used  [58].  Although  there  are  other  accepted  models  for  generating 
random  phase  screens,  including  the  frequently  used  Fourier  Transform  approach  [57],  they 
are  not  addressed  in  this  chapter. 

One  benefit  to  using  the  Zemike  method  is  that  the  Zernike  polynomials  can  be  used 
to  simulate  mirror  aberrations.  Using  Zernike  coefficients  for  multiple  functions  also  helps 
to  reduce  processing  time.  The  Zernike  method  has  a  less  difficult  time  reproducing  the  low 
frequency  effects  of  the  atmosphere,  such  as  tilt,  than  the  Fourier  Transform  method.  The 
Zernike  method  does  not  have  the  same  low  frequency  issues.  The  low  frequency  effects 
have  the  most  power,  and  are  a  dominant  source  of  error  in  most  telescopes.  Additionally, 
high  frequency  aberrations  do  not  have  a  large  impact  on  the  PSF  in  undersampled  systems, 
such  as  the  telescopes  discussed  in  this  dissertation.  High  frequency  changes  will  be  filtered 
by  the  relatively  large  CCD  pixels.  In  this  chapter,  a  total  of  15  Zernike  polynomials  are 
used  to  represent  the  random  aberrations  present  in  the  atmosphere.  Including  a  higher 
number  of  Zernikes  is  not  found  to  have  a  significant  impact  on  the  PSF,  and  as  a  result, 
the  determined  noise  model. 

The  atmosphere  acts  like  a  random  process  to  any  light  that  travels  through  it.  The 
result  is  a  propagation  path  for  the  light  that  has  varying  indices  of  refraction  in  space  and 
time,  n(x,y,z).  This  random  process  causes  distortion  in  the  wavefront,  6a(x,y),  and  can  be 
represented  as  a  sum  of  a  set  of  orthonormal  Zernike  polynomials,  Zj(x,y). 


0a(m,  n )  =  £  bjZj(m,  n ) 


(5.3) 
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where  bj  are  the  weighting  coefficients  for  the  /h  polynomial.  In  this  case,  bj  are  random 
variables  that  represent  the  amount  of  each  type  of  distortion  present  in  the  atmosphere 
that  the  light  is  propagated  through.  The  coefficients  are  zero  mean  random  variables, 
implying  that  there  is  no  tendency  for  any  of  the  coefficients  or  aberrations  to  be  skewed 
in  one  direction  or  another.  The  variance  of  the  Zemikes  due  to  a  Kolmogorov  spectrum, 
as  described  in  chapter  2,  are  derived  in  [62].  These  variances  form  a  covariance  matrix, 
C.  To  utilize  this  covariance  matrix  to  generate  coefficients  with  the  correct  statistics,  a 
Cholesky  factorization  is  performed  on  the  covariance  matrix  C  =  ®®7 .  Defining  n  as  a 
Gaussian  random  vector  with  zero  mean  and  unit  variance,  another  vector  b  =  /7®  based  on 
the  Cholesky  factorization  of  the  covariance  is  formed.  This  vector  b  accurately  represents 
the  Kolmogorov  statistics  present  in  the  atmosphere.  This  can  be  shown  with  the  following 
mean  and  variance  equations: 


E  [£]  =  E  [77®]  =  ®£  [77]  =  0  (5.4) 

E  [bbT]  =  E  [®727?®r]  =  ®£  [7777]  ®r  =  C.  (5.5) 

The  vector  b  is  zero  mean  and  has  covariance  of  C.  b  is  then  used  as  the  weighting 
coefficients  for  the  Zemike  polynomials  described  in  equation  (5.3).  These  atmospheric 
aberrations  are  then  included  in  the  pupil  function,  along  with  the  aperture  and  lens 
aberrations. 

It  is  important  to  note  that  the  phase  screen,  9a(m,  77),  is  one  instance  of  atmosphere 
experienced  by  the  object.  As  mentioned  previously,  the  phase  screen  representing  lens 
aberrations  Q0(m,ri)  is  constant  over  m  and  n.  To  create  images  representing  a  non- 
instantaneous  integration  time,  additional  phase  screens  and  their  resulting  PSFs  need  to 
be  generated.  This  process  is  described  in  section  5.2.1. 
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5.1.2  Detection  Theory. 

One  commonly  used  detection  algorithm  is  a  Log-Likelihood  Ratio  (LLR).  Assuming 
an  Equal-Cost  Equal-Prior  (ECEP)  model,  a  Bayes  risk  criterion  for  a  BHT  can  be 
described  with  the  following  LLR  [25]: 


LLR  =  log 


p(D\Hi)\  H * 
P(D\H0)I 


k  T 

Ho 


(5.6) 


Hi  is  the  hypothesis  that  a  RSO  is  present  in  the  image  and  H0  is  the  hypothesis  that 
there  is  no  RSO  present,  r  is  the  threshold  that  sets  the  acceptable  false  alarm  rate.  D  is 
the  received  data,  which  can  be  either  the  photon  count  in  a  single  pixel  or  across  multiple 
pixels,  D(x,y),  depending  on  if  a  point  detector  or  matched  filter  detector  is  being  used. 
Algorithms  for  the  detection  of  space  objects  have  been  developed  for  both  point  detectors 
[16]  and  matched  filter-based  detectors  [22,  24]  under  different  received  noise  distribution 
assumptions.  The  conditional  probabilities  p(D\Ho)  and  p(D\H\)  are  the  probability  of  the 
received  data  given  that  the  specified  hypothesis  is  true.  The  received  data  is  modeled  with 
the  following  equation: 


D(x,  y)  =  9hn(x,y)  +  B  (5.7) 

D(x,y )  is  the  received  intensity  at  the  CCD  before  photon  counting  noise  is  included. 
The  normalized  PSL  hn(x,y )  is  the  optical  PSL,  hopt(x,y),  normalized  to  sum  to  one.  The 
variable  6  is  the  intensity  of  the  object  being  investigated.  In  the  case  where  an  object  is 
considered  present,  Hi,  6  can  be  any  positive  non-zero  integer.  Under  the  hypothesis  where 
no  object  is  present,  9,  is  zero  and  the  model  reduces  to  a  constant  background  level  B. 
The  conditional  probability  of  data  given  Hi,  p(D\H\ ),  is  the  distribution  of  interest  to  this 
research.  The  distribution  fit  theory  is  described  in  section  5.1.3  and  analyzed  in  section 
5.3.  The  probability  of  the  data  given  H0,  p(D\H0),  is  considered  to  be  Poisson,  due  to  the 
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fact  that  the  dominant  noise  source  present  under  this  hypothesis  is  photon  counting  noise 
[17].  The  Poisson  conditional  probability  can  be  expressed  with  the  following  equation: 


BDb  e~B 

p(DB\H0)  =  — — .  (5.8) 

uB'- 

Db  is  the  number  of  background  photons  received,  and  B  is  the  expected  number  of 
background  photons. 

5.1.3  Received  Noise  Distributions. 

As  mentioned  previously,  detection  algorithms  in  electro-optical  telescopes  focus 
primarily  on  two  types  of  PDF  for  noise  statistics,  Gaussian  [24,  56]  and  Poisson  [16,  35]. 
To  justify  these  distributions,  limiting  assumptions  are  made.  The  assumptions  made 
can  degrade  or  limit  the  performance  of  the  detection  algorithm  and  the  capabilities  of 
the  telescope  system.  Often  these  noise  assumptions  demonstrate  a  trade-off  between 
computational  complexity  and  increased  performance.  For  example,  one  benefit  of  a 
Gaussian  noise  assumption  is  that  it  allows  for  easier  computation  of  the  SNR.  This  may 
potentially  result  in  reduced  performance,  if  the  received  data  is  not  truly  Gaussian. 

The  model  proposed  in  this  chapter  combines  two  important  effects  into  a  mixed 
distribution  to  attempt  to  achieve  a  closer  match  to  the  true  measured  noise  distribution.  In 
previous  research  efforts  [36,  56],  a  Gamma  distribution  is  investigated,  both  theoretically 
and  analytically,  to  describe  the  statistics  of  intensity  fluctuations  in  the  atmosphere.  There 
have  been  other  proposed  distributions,  including  Log-normal  and  Beta  [56],  but  these  are 
not  explored  in  this  chapter. 

The  derivation  of  the  mixed  PDF  models  begins  with  the  assumption  that  the 
underlying  intensity  fluctuations  are  Gamma  distributed.  Once  the  perturbed  field  from  the 
space  object  reaches  the  CCD  and  is  captured,  the  photons  experience  Poisson  counting 
noise.  Mathematically  this  combination  of  distributions  is  achieved  by  combining  the 
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continuous  Gamma  function  with  the  discrete  photon  counting  noise  present,  leading  to 
a  negative  binomial  distribution  for  p(D\H\)  [72]: 


P(Ds\Hi) 


y(m  +  ds) 


1  +  — 

M 


-M 


1  +  — ) 
s  1 


-Ds 


(5.9) 


r(M)r(Ds) 

r  is  the  Gamma  function  and  M  is  the  number  of  degrees  of  freedom.  M  represents  the 
amount  of  scintillation  present  in  the  received  data.  A  value  of  one  for  M  gives  the  largest 
scintillation  effect.  Alternatively,  larger  M  values  decrease  the  amount  of  scintillation 
present  in  the  received  data.  S  is  the  intensity  of  the  RSO  of  interest,  and  D,  is  the  number 
of  photons  received  from  the  target  object. 

The  scintillated  received  data  described  to  this  point  does  not  include  the  total  number 
of  photons  captured.  In  addition  to  photons  for  the  object  being  observed,  background 
photons  from  ambient  light  sources  will  also  be  captured  in  the  data.  The  total  number 
of  photons  received  by  the  CCD,  D,  is  a  combination  of  the  number  of  photons  from  the 
object  and  the  number  of  photons  from  the  background:  D  =  D,  +  DB.  Assuming  that 
background  and  scintillated  light  from  RSO  are  statistically  independent,  the  two  sources 
can  be  combined  into  a  joint  Probability  Mass  Function  (PMF).  To  achieve  a  marginal 
PMF,  the  joint  PMF  is  summed  over  the  number  of  photons  received  from  the  object  Ds 
under  the  H\  hypothesis: 


P(D\H{)  = 


e~BBD  /  |  S  \~M  y  Y(M  +  D,)  /  |  Af\“D'  B°s 

T(M)  l1  +  M)  ^  Y(DS)  V  +  j)  (D-D,)! 


(5.10) 


There  are  two  main  challenges  presented  by  the  mixed  PMF  shown  in  5.10.  The 
first  is  that  the  conditional  probability  given  H\  is  not  easily  reduced  to  a  simple  PDF. 
Additionally,  when  combined  to  get  a  LLR,  the  new  method  does  not  reduce  to  a  sufficient 
statistic  similar  to  Gaussian  and  Poisson.  Instead  of  a  simple  sufficient  statistic  such  as 
SNR  found  through  a  correlation  of  the  data  and  expected  PSF,  the  probability  needs  to 
be  calculated  directly.  Secondly,  the  distribution  relies  on  two  parameters  that  are  not 
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immediately  apparent  from  the  received  data:  the  coherence  parameter  M  and  the  true 
object  intensity  S .  These  parameters  need  to  be  estimated  from  the  data.  Methods  of 
determining  these  parameters  are  described  in  section  5.3.1. 

5.2  Simulation  Description 

In  this  section,  a  simulation  is  developed  to  generate  realistic  data  that  would  model  the 
received  data  from  a  telescope  observing  a  RSO.  The  model  accounts  for  both  scintillation 
from  the  atmosphere  and  photon  counting  noise  effects.  The  simulated  data  serves  two 
purposes  in  this  research.  The  first  is  to  analyze  the  intensity  fluctuations  present  in  the 
PSF,  and  determine  which  PMF  is  the  best  fit.  After  the  most  appropriate  distribution  fit 
has  been  determined,  the  RSO  data  is  processed  with  the  detection  algorithm  developed 
in  the  previous  section  to  evaluate  the  performance  against  traditionally  used  detection 
algorithms. 

5.2.1  Evolving  Phase  Screen  Over  Time. 

The  weighting  vector  coefficients,  b.  developed  in  the  previous  section  accurately 
represent  the  statistics  of  the  atmosphere  at  one  instant  in  time.  To  simulate  any  non- 
instantaneous  and  measurable  period  of  exposure  time,  the  behavior  of  the  atmosphere 
over  time  is  essential.  It  has  been  shown  that  unless  there  are  very  strong  winds,  or  a  long 
time  between  observations,  the  statistics  of  the  atmosphere  are  correlated  over  time  [58]. 
This  implies  that  over  a  short  period  of  time,  the  optical  system  is  likely  to  see  correlated 
phase  errors,  6a(m,  n )  across  the  pupil.  It  is  important  to  account  for  this  effect  in  order  to 
properly  simulate  the  behavior  over  time.  To  account  for  the  correlated  atmosphere,  each 
of  the  j  Zernikes  modeled  have  a  corresponding  correlation  term,  Rj.  Taylor’s  frozen  flow 
hypothesis  is  used  to  create  this  model  of  correlation.  In  Taylor’s  Frozen  Flow,  a  constant 
wind  blows  a  frozen  atmosphere  across  the  pupil  [48].  This  correlation  for  the  jth  Zemike 
is  a  function  of  the  separation  between  two  points  in  time  along  with  wind  velocities  and 
seeing  parameter  ra. 
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Rjm  =  E[bj(h)bj(t2)\ 


(5.11) 


Examining  the  expected  value  of  the  coefficient  at  these  two  times  shows  how 
correlated  each  new  Zemike  at  time  t2,  bj(t2),  will  be  to  the  previous  at  time  t\,  bj{t{). 
The  correlation  for  a  specific  time  separation,  A t  =  t2  —  t\,  is  found  numerically  through 
the  structure  function  of  the  atmosphere.  As  screens  evolve  over  time,  the  weighting 
coefficients  are  no  longer  zero  mean,  and  have  a  variance  of  cra..  A  new  mean  is  calculated 
based  on  the  previous  values  of  each  coefficient  bj,  and  the  correlation  R /(At). 


A/  =  bj 


R,m 


CT  ■ 


(5.12) 


The  variance  of  the  next  instance  of  Zemike  coefficients  is  also  a  conditional  variance. 
The  variance  is  based  on  the  unconditioned  variance  of  the  Zernike  and  the  correlation  for 
the  specified  time  separation  of  At. 


= 


W)2  - 


cr2 
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(5.13) 


A  new  vector  of  coefficients,  c;,  is  formed  from  the  conditional  mean  and  variance 
in  equations  (5.12)  and  (5.13).  These  coefficients  are  also  Gaussian  distributed  random 
variables,  but  with  a  mean  fij  and  variance  &/. 


Cj~N(m,d-j)  (5.14) 

The  variance  in  the  fluctuations  is  related  to  the  ratio  between  the  diameter  of  the 
aperture  of  the  telescope,  d  and  Fried’s  Atmosphere  coherence  diameter,  r0  [17].  This 
ratio,  d/r0,  determines  the  amount  of  each  type  of  phase  error  present. 

5.2.2  Final  System  Model  and  Implementation. 

To  finalize  the  model  for  the  simulation,  the  telescope  and  random  phase  screen 
simulations  are  combined.  There  are  different  approaches  to  simulating  propagation 
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through  the  atmosphere  and  placement  of  the  phase  screen.  Figure  5.1  demonstrates  two 
distinct  methods  of  simulation. 


Figure  5.1:  Diagram  of  two  potential  different  atmosphere  propagations.  In  (A)  a  multiple 
step  propagation  is  demonstrated,  in  (B)  a  single  phase  screen  directly  against  the  pupil  is 
demonstrated. 


In  (A),  the  field  is  propagated  through  multiple  phase  screens  placed  along  the 
propagation  path.  Each  screen  may  contain  different  statistics  that  represent  the  section 
of  atmosphere  around  that  location.  This  method  allows  for  more  complicated  interactions 
between  the  field  after  each  new  simulated  section  of  atmosphere.  In  (B),  the  phase  screen, 
6a(x,y),  is  simulated  to  be  a  thin  phase  screen  placed  directly  against  the  pupil  of  the 
telescope.  Generally  this  method  may  not  accurately  represent  the  atmospheric  effects, 
unless  an  isoplanatic  assumption  is  valid.  The  isoplanatic  assumption  implies  that  light 
from  any  part  of  an  object  of  interest  experiences  the  same  section  of  atmosphere.  To 
ensure  this  is  a  valid  assumption,  a  model  for  the  atmosphere  is  needed. 
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One  popular  model  of  refractive-index  structure  constant  C\  is  Hufnagal- Valley 5/7. 
This  structure  constant  is  a  measure  of  the  strength  of  fluctuations  in  the  refractive  index. 
Under  normal  vertical  viewing  conditions,  it  gives  an  isoplanatic  angle  of  7/rrad  [36].  The 
isoplanatic  angle  is  an  angle  over  which  the  isoplanatic  assumption  can  be  considered  valid. 
Most  SDA  search  telescopes  have  a  large  field  of  view,  much  much  larger  than  7/rrad,  so 
an  isoplanatic  assumption  is  not  valid  when  simulating  an  entire  field  of  view.  However, 
since  only  a  single  point  object  is  being  simulated  and  investigated  to  generate  a  PSF, 
an  isoplanatic  assumption  can  be  accurately  made.  For  example,  a  150m  object  in  GEO 
would  still  only  represent  an  angle  of  4.2  //rad,  which  is  less  than  the  isoplanatic  angle  in 
the  Hufnagal- Valley  model.  Now  that  there  is  a  model  for  the  phase  screens,  9a(m,  n),  these 
simulated  phase  errors  created  from  equation  (5.14)  are  added  to  the  pupil  function  from 
equation  (5.2)  to  give  an  updated  pupil  function,  Pa(m,  n),  that  is  used  to  generate  the  PSF. 

Pa(m,  n )  =  Aim,  n)  exp  [j90(m,  n)\  exp  \jdjm,  77)]  (5.15) 

This  model  now  accounts  for  the  first  phenomenon  being  considered,  the  scintillation 
due  to  the  atmosphere.  The  next  step  is  to  add  in  the  background  light  as  shown  in  equation 
(5.7).  Finally,  photon  counting  Poisson  noise  is  added  on  top  of  the  expected  image  D(x,y ) 
to  give  the  final  received  data. 

The  optical  model,  random  phase  screens,  and  data  analysis  are  processed  in 
MATLAB.  All  phase  screens  are  generated  with  a  separation  of  At  =  1ms.  To  simulate 
a  desired  exposure  time,  1/A t  correlated  phase  screens  are  used  and  each  is  put  through 
the  combined  model.  Since  a  non-instantaneous  exposure  time  is  desired,  it  is  necessary 
to  evolve  the  Zemike  polynomials  coefficients  to  generate  each  subsequent  phase  screen, 
as  described  in  section  5.2.1.  In  each  new  Monte  Carlo  simulation,  an  independent  phase 
screen  is  generated.  Next,  the  remaining  (1/At  -  1)  phase  screens  are  generated  based  on 
the  statistics  of  the  first  screen  and  the  correlation  of  Zernike  polynomials.  Each  of  these 
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screens  are  then  placed  in  the  optical  model  to  create  PSFs.  The  resulting  PSF  images  are 
averaged  and  normalized  to  create  a  single  simulated  PSF.  This  PSF  is  the  result  of  the 
desired  exposure  time,  optical  system,  and  atmospheric  conditions. 

5.3  Results  and  Analysis 

To  demonstrate  that  the  received  data  more  accurately  matches  the  newly  developed 
mixed  PDF  than  traditional  assumption,  the  mixed  PDF  is  compared  against  both  Gaussian 
and  Poisson  distributions.  An  error  metric  is  used  to  determine  the  best  fit  between  the 
simulated  received  data  and  the  potential  distributions  for  various  atmospheric  conditions. 

5.3.1  Distribution  Matching. 

The  first  step  in  analyzing  the  distribution  fits  is  to  create  a  histogram  that  represents 
the  PDF  of  the  received  data.  This  PDF  can  then  be  compared  against  the  theoretical 
PDF  of  the  potential  distributions.  The  center  pixel  intensity  values  in  the  PSF  are  noted 
across  1000  independent  Monte  Carlo  iterations  and  then  normalized.  An  important  aspect 
in  fitting  the  distributions  is  determining  the  correct  parameters  for  each  PDF  based  on 
the  received  intensities.  The  Gaussian  is  found  by  investigating  the  mean  and  standard 
deviation  of  the  received  data.  The  Poisson  is  described  completely  by  the  mean.  The 
mixed  PDF  depends  on  four  parameters:  M,  S ,  B,  and  I).  The  value  of  D  is  the  simulated 
data  of  the  pixel  being  investigated.  The  value  of  M  is  found  by  providing  the  true  values 
for  S  and  B ,  then  selecting  the  M  that  provides  the  best  fit.  The  value  of  M  is  a  function  of 
atmospheric  conditions  and  the  integration  time  Ts.  Once  the  fit  for  M  is  found,  the  values 
for  S  and  B  are  estimated  from  the  data. 

Each  PDF  is  evaluated  against  the  received  PDF  at  the  histogram  bin  locations.  The 
difference  between  PDFs  is  considered  the  error  between  simulated  data  and  the  candidate 
PDF.  Each  vector  of  errors  is  then  combined  into  a  total  RMSE,  ERMS: 
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(5.16) 


E RMS  - 

where  N/,  is  the  number  of  buckets,  or  places  in  the  histogram  that  the  PDFs  are  being 
compared  against.  /  is  the  bin  locations,  p(l )  is  the  value  of  the  PDF  being  investigated  at 
bin  /,  and  D(l)  is  the  simulated  data  value  evaluated  at  the  histogram  bin  /.  Figure  5.2  shows 
the  £rms  values  for  the  three  candidate  PDFs  considered  versus  the  seeing  parameter,  ra. 
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Figure  5.2:  Error  between  candidate  PDFs  and  received  intensity  fluctuations.  £rms  is 
calculated  for  a  variety  of  seeing  values,  r0 


The  error  between  the  mixed  distribution  and  the  simulated  data  is  smallest  for  r0 
values  less  than  6  cm.  After  6cm,  the  Gaussian  PDF  more  closely  matches  the  received 
data.  Physically  this  means  that  during  conditions  of  increased  turbulence,  or  lower  rQ,  the 
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temporal  scintillation  will  be  more  prominent.  There  are  several  locations  where  survey 
telescopes  are  located  that  routinely  experience  r0  values  of  less  than  6cm,  indicating  that  it 
is  not  unrealistic  to  expect  these  conditions.  It  is  found  that  the  distribution  match  is  largely 
a  factor  of  the  seeing  parameter,  ra,  used.  Additional  factors  including  the  integration 
time,  wind  speed,  and  number  of  Zemikes  used  demonstrate  a  negligible  impact  on  the 
distribution  fit.  Although  the  Poisson  does  not  have  the  lowest  error  for  the  data  analyzed, 
it  can  not  be  universally  stated  that  Gaussian  is  better  than  Poisson.  This  data  set  was 
generated  making  assumption  about  the  noise  sources  present.  Not  all  possible  noise 
sources  are  considered. 

These  findings  indicate  that  for  certain  conditions,  the  mixed  PDF  results  in  the  closest 
distribution  match  to  the  simulated  PSF  fluctuations.  An  important  note  is  that  the  analysis 
performed  is  based  on  the  statistics  of  a  single  pixel,  the  center  pixel  of  the  PSF.  To 
generalize  this  analysis  to  include  the  entire  PSF,  it  is  assumed  that  each  pixel  in  the 
PSF  is  independent.  Moving  from  a  point  detection  strategy  to  a  matched  filter,  or  entire 
PSF  strategy,  allows  for  more  advanced  detection  algorithms.  In  reality,  there  may  be  a 
correlation  in  the  temporal  scintillation  of  adjoining  pixels  that  propagate  through  similar 
portions  of  the  atmosphere.  This  effect  is  not  analyzed  in  this  research  effort.  In  the  next 
section,  the  mixed  distribution  is  implemented  in  a  LLR  test  to  determine  the  performance 
of  the  newly  developed  algorithm.  The  algorithms  are  compared  on  their  ability  to  find 
potential  RSOs. 

5.3.2  Receiver  Operating  Characteristic  Curves. 

A  Receiver  Operating  Characteristic  (ROC)  curve  is  a  method  of  comparing  the 
detection  performance  across  the  entire  range  of  false  alarm  values.  It  is  generated  by 
sliding  the  threshold  r  from  equation  (5.6),  and  noting  Pc\  and  Pf  value  pairs  at  that 
threshold.  To  demonstrate  the  maximum  performance  achievable  with  this  detection 
algorithm,  simulated  data  is  created  with  the  mixed  PDF  developed  in  section  5.1.3. 
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Figure  5.3  shows  two  instances  of  simulated  PSFs  created  with  (A)  a  long  exposure 
atmosphere,  and  (B)  an  atmosphere  scintillated  with  the  mixed  negative  binomial  and 
Poisson  distribution. 


Figure  5.3:  A  demonstration  of  a  simulated  received  PSF  under  two  different  conditions. 
In  (A)  a  long  exposure  atmosphere  model  is  used,  and  in  (B)  a  PSF  generated  with  the 
negative  binomial  and  Poisson  model  or  the  ideal  mixed  PDF. 


Using  a  Monte  Carlo  simulation  of  the  PSFs,  the  detection  performance  of  the  two 
algorithms  that  most  closely  match  the  PSF  fluctuations  are  investigated.  These  two 
distributions  are  the  mixed  and  Gaussian  PDFs.  To  produce  a  ROC  curve,  two  types  of 
data  are  needed.  These  two  data  sets  are  images  with  and  without  a  RSO.  Performing 
the  detection  algorithm  on  data  with  an  object  present  gives  the  probability  of  detection, 
while  performing  the  detection  algorithm  on  data  with  no  object  present  determines  the 
probability  of  false  alarm.  Images  without  an  RSO  are  generated  by  creating  a  flat 
background  with  a  value  of  B,  and  adding  Poisson  photon  counting  noise. 

One  method  of  determining  Pd  and  Pf  is  to  compute  the  LLR  and  compare  it  against 
a  threshold  and  assign  either  a  one  for  a  detection  or  a  zero  for  no  detection.  This  can  be 
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done  over  a  large  number  of  independent  Monte  Carlo  iterations.  The  ratio  of  detection 
or  false  alarms  to  the  number  of  Monte  Carlo  iterations  gives  a  value  for  Pd  and  Pf.  This 
direct  computational  method  requires  significant  data  processing  at  each  threshold  value. 
In  addition,  the  ROC  curve  is  limited  in  false  alarm  range  by  the  number  of  iterations 
performed.  To  achieve  false  alarm  rates  of  10“5,  on  average  105  iterations  would  be  needed 
to  generate  a  single  false  alarm. 

Instead,  a  Monte  Carlo  simulation  is  performed  with  a  large  number  of  iterations, 
N  =  1000,  and  the  LLR  statistics  are  analyzed.  By  investigating  the  statistics  of  the  LLR, 
values  for  Pd  and  Pj  can  be  calculated  with  a  CDF.  Figure  5.4  shows  N  =  100  instances  of 
LLR  generated  under  two  conditions  (A)  where  the  Gaussian  detector  is  implemented  and 
(B)  where  the  mixed  PDF  detector  is  implemented. 

(A)  (B) 


Figure  5.4:  LLR  observations  over  100  Monte  Carlo  iterations.  The  LLR  is  generated  for 
both  the  (A)  Gaussian  algorithm  and  the  (B)  mixed  PDF  algorithm. 


The  LLR  values  in  the  case  where  no  RSO  is  present,  shown  as  the  dashed  line  in 
Figure  5.4,  for  both  Gaussian  and  mixed  PDF  algorithms,  closely  follows  a  Gaussian 
distribution.  Due  to  this  fact,  the  LLR  statistics  can  be  summarized  by  their  mean  and 
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standard  deviation.  The  LLR  values  with  an  RSO  present  utilizing  the  Gaussian  detection 
algorithm  is  also  Gaussian.  The  exceptions  are  the  LLR  values  with  an  object  present  in  the 
mixed  PDF  algorithm.  Analyzing  the  data  shows  that  it  more  closely  follows  an  exponential 
distribution.  This  can  be  seen  in  Figure  5.4  (B),  as  the  data  tends  to  skew  upwards  of  the 
mean  instead  of  symmetrically  around  the  mean. 

In  Figure  5.4,  Pd  is  found  by  computing  the  CDF  of  the  solid  line  using  the  statistics  of 
the  LLR  and  a  sliding  threshold  value  r.  Alternatively,  Pf  is  found  by  computing  the  CDF 
using  the  statistics  of  the  dashed  line  and  the  threshold  r.  To  ensure  that  the  entire  range 
of  Pd  and  Pf  values  are  investigated,  it  is  necessary  to  use  a  threshold  that  spans  from  the 
most  negative  to  the  most  positive  potential  LLR  values.  Using  this  method,  ROC  curves 
are  found  for  both  algorithms,  and  are  shown  in  Figure  5.5. 

For  high  Pj  values,  the  algorithms  perform  similarly,  but  the  mixed  PDF  algorithm 
shows  Pd  improvement  as  the  acceptable  Pf  rate  is  decreased.  The  new  algorithm 
demonstrates  detection  improvements  of  up  to  35  percent  for  sufficiently  low  false  alarm 
rates.  The  mixed  PDF  algorithm  appears  to  begin  at  approximately  Pf  =  0.01,  but  actually 
has  a  Pd  value  of  one  and  is  difficult  to  distinguish  from  the  plot  border.  Low  Pf  rates  are 
typically  where  SDA  systems  operate,  and  as  a  result  this  algorithm  can  potentially  have  a 
benefit  over  the  Gaussian. 

For  telescope  systems  with  CCDs  containing  a  large  number  of  pixels,  a  low  Pf  is 
desirable.  For  example,  a  512x512  array  utilizing  a  false  alarm  rate  of  10-5  would  still 
average  approximately  2.6  false  alarms  per  image.  Depending  on  the  program  objectives 
and  process  timing  requirements,  this  may  or  may  not  be  an  acceptable  number  of 
incorrectly  identified  RSOs.  The  performance  achieved  in  this  ROC  curve  is  the  theoretical 
maximum  improvement.  Factors  including  the  independence  of  adjacent  pixels  and  how 
closely  the  atmospheric  conditions  match  the  simulated  parameters  will  determine  the 
realizable  detection  increase  of  the  new  algorithm,  including  how  many  more  RSOs  can 
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Figure  5.5:  ROC  curves  comparing  the  new  mixed  PDF  (solid  line)  against  the  traditional 
Gaussian  detection  algorithm  (dashed  line).  These  curves  were  generated  with  simulated 
RSO  and  background  data. 


be  found.  The  better  the  seeing  parameter  ra  experienced,  the  less  effective  the  new  mixed 
PDF  algorithm  would  perform.  The  equal  performance  point  occurs  at  approximately  6cm, 
and  the  Gaussian  algorithm  has  a  larger  Pd  for  larger  r0  values. 

Additionally,  factors  including  spatial  aliasing  and  integration  time  have  a  large  impact 
on  the  effectiveness  of  this  mixed  PDF  model.  Spatial  averaging  due  to  aliasing  can  hide 
the  intensity  fluctuations  that  give  this  new  algorithm  its  benefits.  In  the  case  of  long 
integration  times,  the  phase  screen  model  is  less  accurate,  and  a  long  exposure  atmosphere 
more  accurately  represents  the  received  data. 
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5.4  Conclusions 


This  research  demonstrates  that  under  certain  atmospheric  conditions,  a  mixed  PDF 
model  combining  negative  binomial  scintillation  and  Poisson  background  light  can  more 
accurately  represent  the  received  noise  statistics  of  the  atmosphere  present  in  an  observed 
PSF  than  a  Gaussian  assumption.  For  seeing  values  of  less  than  6cm,  this  new  PDF  is  found 
to  more  accurately  match  received  intensity  fluctuations.  Using  this  model  under  ideal 
assumptions,  the  mixed  PDF  model  shows  improvements  of  up  to  35  percent  detection  for 
realistic  Pf  values  typically  used  in  SDA  telescopes.  Accounting  for  the  scintillation  that 
is  present  in  the  RSO  data,  but  not  the  background  data  helps  distinguish  objects  from  the 
background.  In  this  way,  the  additional  atmospheric  noise  improves  detection  performance 
in  this  model. 

There  are  a  few  areas  related  to  the  research  presented  in  this  chapter  where  additional 
future  work  could  be  performed.  This  work  could  be  expanded  in  the  future  by  applying  the 
algorithm  to  collected  telescope  data.  Only  simulated  data  is  investigated  in  this  chapter, 
and  to  enable  future  implementation  it  is  important  to  test  the  algorithms  on  collected 
telescope  data.  There  are  other  methods  of  hypothesis  testing  besides  LLR.  Using  an 
alternative  detection  algorithm  along  with  the  noise  model  described  in  this  chapter  may 
show  additional  benefits  over  a  LLR. 

This  research  did  not  investigate  the  optimal  methods  of  estimating  object  signal 
strength  S  and  the  background  B.  There  may  be  benefit  to  optimizing  the  estimation 
techniques  for  these  parameters.  In  addition,  it  is  assumed  that  the  pixels  in  the  data  are 
statistically  independent  as  well  as  the  background  and  photons  from  the  object  of  interest. 
The  model  could  be  expanded  further  by  investigating  the  impacts  of  these  assumptions 
being  incorrect. 

The  next  chapter  discusses  the  conclusions  and  results  from  throughout  this 
dissertation,  as  well  as  combining  areas  of  future  research.  The  contributions  of  this 
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research  effort  as  well  as  a  broad  overview  of  how  research  fits  into  the  field  is  also 
discussed. 
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VI.  Summary  and  Future  Work 


This  chapter  will  summarize  results  and  contributions  of  this  dissertation  as  well  as 
present  potential  topics  for  additional  future  work. 

6.1  Work  Completed 

In  chapter  1,  a  research  goal  and  three  associated  research  questions  are  presented. 
Throughout  this  dissertation  and  research,  each  of  these  questions  is  addressed  and 
answered.  In  the  following  list,  a  summary  of  the  research  results  that  answer  each  question 
is  included. 

1.  Will  a  new  realistic  cost  function  improve  the  detection  performance  of  a  Bayes 
Criterion  MHT? 

In  chapter  3,  it  is  shown  that  a  new  cost  function  in  a  MHT  improves  detection 
performance  over  traditional  point  detection  approaches,  BHT  algorithms  like  Source 
Extractor,  and  ECEP  MHT  algorithms.  It  is  demonstrated  that  the  probability  of 
detection  is  increased  by  using  an  UCEP  MHT  detection  algorithm  when  compared 
to  an  ECEP  algorithm.  Pd  gains  of  up  to  80  percent  are  observed  for  the  same  Pf 
rate.  The  number  of  hypotheses  used  within  the  MHT  is  also  investigated.  It  is  found 
that  there  is  a  relationship  between  the  number  of  hypotheses  used  and  detection 
performance.  The  results  show  that  a  MHT  clearly  detects  more  objects  than  a  BHT, 
and  that  using  6  hypotheses  slightly  outperforms  the  algorithm  when  10  hypotheses 
are  used,  with  the  additional  benefit  of  reduced  computational  time. 

2.  Do  the  assignments  of  a  priori  probabilities  in  a  MHT  improve  the  detection 
performance  ? 

Chapter  4  addresses  this  research  question.  It  is  shown  that  increased  performance 
is  achievable,  in  both  simulated  and  collected  SST  data.  When  tested  with  simulated 
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data,  ROC  curve  analysis  shows  that  the  ECUP  algorithm  has  a  greater  Pd  than  the 
ECEP  algorithm  for  a  large  range  of  Pf  values.  At  a  relevant  false  alarm  level  of 
Pf  =  10“10,  ECUP  shows  an  average  improvement  of  approximately  5  percent  when 
tested  over  100  Monte  Carlo  iterations. 

When  tested  against  collected  SST  data,  this  research  demonstrates  that  for  certain 
applications  it  may  be  beneficial  to  use  an  ECUP  MHT  algorithm.  It  is  shown  with 
high  statistical  confidence  that  the  ECUP  algorithm  on  average  has  more  detections 
than  an  ECEP  algorithm  when  performed  on  the  SST  data.  For  the  small  subset 
of  frames  analyzed,  approximately  19-26  additional  detections  are  found,  which 
translates  to  458-627  additional  detections  on  average  for  full  frame  data.  These 
objects  are  most  likely  threshold  objects  that  are  difficult  for  currently  proposed 
detection  algorithms  to  detect.  Finding  these  threshold  objects  will  provide  additional 
information  for  tracking  and  characterization,  in  turn  increasing  SDA  capabilities 
with  the  optical  system. 

3.  Can  using  a  more  realistic  noise  model  for  detection  algorithms  increase  the  ability 
to  detect  space  objects? 

In  chapter  5  it  is  shown  that  under  certain  atmospheric  conditions,  a  mixed  PDF 
model  combining  negative  binomial  scintillation  and  Poisson  background  light  can 
more  accurately  represent  the  received  noise  statistics  of  the  atmosphere  present  in  an 
observed  PSF  than  a  Gaussian  assumption.  For  seeing  values  of  less  than  6cm,  this 
new  PDF  is  found  to  more  accurately  match  received  intensity  fluctuations.  Using 
this  model  under  ideal  assumptions,  the  mixed  PDF  model  shows  improvements  of 
up  to  35  percent  detection  for  realistic  Pf  values  typically  used  in  SDA  telescopes. 
Accounting  for  the  scintillation  that  is  present  in  the  RSO  data,  but  not  the 
background  data,  helps  distinguish  objects  from  the  background.  In  this  way,  the 
additional  atmospheric  noise  improves  detection  performance  in  this  model. 
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6.2  Algorithm  Comparison 

Throughout  this  dissertation,  several  algorithms  are  presented  that  addresses  the  issues 
of  spatial  undersampling  in  CCD  arrays,  as  well  as  the  assumed  noise  distribution.  In  this 
section,  the  performances  of  the  algorithms  developed  are  compared  and  contrasted. 

In  a  MHT,  there  are  two  key  variables:  the  cost  and  prior  probabilities.  Table  6.1 
shows  how  the  chapters  of  this  dissertation,  as  well  as  a  previous  research  effort,  relate  to 
choices  for  costs  and  priors. 


Table  6.1:  Description  of  MHT  algorithms  including  the  assignment  of  costs  and  prior 
probabilities  and  where  they  are  developed. 


Equal- Cost 

Unequal-Cost 

Equal-Prior 

Zingarelli  et.  al.  [24] 

Chapter  3:  UCEP 

Unequal-Prior 

Chapter  4:  ECUP 

N/A 

Table  6.1  shows  research  into  four  different  assumptions  that  can  be  used  in  a 
hypothesis  test.  In  this  case,  the  research  relates  to  improvements  in  SDA  detection 
techniques.  Figure  6.1  is  a  diagram  of  the  general  relationship  between  computational 
complexity  and  Pci  for  SDA  algorithms  from  other  sources,  as  well  as  newly  developed 
algorithms  proposed  in  this  dissertation. 

Moving  along  the  y-axis  in  figure  6.1  from  bottom  to  the  top  increases  the  probability 
of  detection.  The  algorithms  above  the  x-axis  all  utilize  multiple  pixels  to  decide  if  an 
object  is  present.  The  point  detector  on  the  other  hand  uses  only  a  single  pixel  to  make 
decisions.  Moving  along  the  x-axis  from  left  to  right  increases  the  complexity  of  the 
algorithm.  Algorithms  to  the  right  of  the  y-axis  use  a  MHT,  while  the  matched  filter  and 
point  detection  techniques  use  a  BHT. 
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Figure  6.1:  Pc\  versus  computational  complexity  diagram  for  several  SDA  algorithms. 
Processing  time  for  a  full  frame  of  SST  data  is  also  included. 


The  ECUP  algorithm  developed  in  Chapter  4  provides  an  improvement  in  Pd,  and 
is  slightly  more  complex  due  to  the  calculation  of  the  weighting  term  W*.  The  UCEP 
algorithm  detects  even  more  objects,  but  is  more  computationally  complex. 

6.3  Future  Research  Items 

There  are  several  opportunities  for  continued  research  based  on  the  work  presented 
in  this  dissertation.  In  each  chapter  ideas  are  included  for  potential  future  research.  This 
section  compiles  and  describes  these  future  research  items. 
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Chapter  3  future  work: 

This  research  focuses  on  received  data  that  is  assumed  to  be  Gaussian.  Through  this 
research,  a  method  for  implementing  a  different  model  for  received  data  is  developed 
that  can  be  used  on  any  potential  distribution  of  received  data.  The  limiting  cases  of 
costs  are  considered  in  this  chapter.  There  are  many  additional  cost  structures  in  between 
the  two  presented  here.  For  specific  applications,  these  cases  may  give  better  detection 
performance.  The  drawback  is  that  they  will  not  reduce  as  well  as  the  equal-  or  unequal- 
cost  algorithms  described  in  this  dissertation.  Additionally,  the  unequal  cost  method  here  is 
more  computationally  complex  than  the  currently  implemented  algorithms.  Simplification 
and  efficient  implementation  are  not  a  focus  of  this  research.  To  use  this  algorithm 
operationally,  these  factors  will  need  to  be  investigated  further. 

Chapter  4  future  work: 

The  ECUP  algorithm  developed  assumed  Gaussian  received  noise  in  this  chapter.  In  the 
future,  removing  the  Gaussian  assumption  may  detect  additional  space  objects.  Additional 
methods  can  be  combined  to  potentially  improve  detection  performance  even  further.  For 
example,  using  a  Unequal-Cost  Unequal-Prior  (UCUP)  algorithm  may  improve  detection 
performance  when  compared  to  an  ECUP  algorithm. 

In  the  ECEP  MHTs,  the  threshold  is  considered  constant  for  each  hypothesis.  As  a 
result,  the  modeled  intensity  is  effectively  different  for  each  hypothesis,  dependent  on  the 
distribution  of  the  PSF.  This  assumption  is  changed  in  this  paper  for  the  ECUP  algorithm. 
This  assumption  can  be  removed  in  a  future  ECEP  test,  and  may  provide  interesting 
research  opportunities  and  advancements. 

Another  consideration  is  how  including  unequal-priors  may  affect  the  ability  to 
translate  the  detected  objects  into  accurate  tracking  and  orbit  determination.  The  ECUP 


124 


algorithm  takes  into  account  the  probability  of  the  alternate  hypotheses,  each  of  which 
correspond  to  a  sub-pixel  position.  The  sub-pixel  information  gathered  from  performing  a 
ECUP  test  may  lead  to  improved  sub-pixel  position,  and  in  turn,  more  accurate  tracking  of 
detected  RSOs. 

Chapter  5  future  work: 

There  are  a  few  areas  related  to  the  research  presented  in  this  chapter  where  additional 
future  work  could  be  performed.  This  work  could  be  expanded  in  the  future  by  applying 
the  algorithm  to  collected  telescope  data.  Only  simulated  data  was  investigated  in  this 
chapter.  To  enable  future  implementation,  it  is  important  to  test  the  algorithms  on  collected 
telescope  data.  There  are  other  methods  of  hypothesis  testing  besides  LLR.  Using  an 
alternative  detection  algorithm  along  with  the  noise  model  described  in  this  chapter  may 
yield  additional  benefits  over  a  LLR. 

This  research  did  not  investigate  the  optimal  methods  of  estimating  object  signal 
strength  S  and  the  background  B.  There  may  be  benefits  to  optimizing  the  estimation 
techniques  for  these  parameters.  In  addition,  it  is  assumed  that  the  pixels  in  the  data  are 
statistically  independent,  as  well  as  the  background  and  photons  from  the  object  of  interest. 
The  model  could  be  expanded  further  by  investigating  the  impacts  of  these  assumptions 
being  incorrect. 

One  final  area  that  may  be  beneficial  to  address  in  the  future  is  an  Unequal-Cost 
Unequal-Prior  (UCUP)  algorithm.  This  algorithm  is  not  discussed  in  this  dissertation,  but 
could  improve  detection  performance  above  what  was  demonstrated  with  UCEP.  It  would 
be  the  most  computationally  complex,  and  would  not  provide  sub-pixel  information,  as  is 
the  case  with  the  UCEP  algorithm. 
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6.4  Publications 


To  ensure  this  work  is  novel  and  well  received  by  the  community,  several  portions 
of  this  dissertation  have  been  submitted  for  conference  and  journal  papers.  In  total,  two 
journal  articles  and  two  conference  papers  have  been  written  and  submitted. 

Journals: 

•  Based  on  Chapter  3: 

T.  Hardy,  S.  Cain,  J.  Jeon,  and  T.  Blake.  “ Improving  space  domain  awareness 
through  unequal-cost  multiple  hypothesis  testing  in  the  space  surveillance  telescoped 
Applied  Optics  54.17  (2015):  5481-5494. 

•  Based  on  Chapter  4: 

T.  Hardy,  S.  Cain,  and  T.  Blake.  “ Unequal  A  Priori  Probability  Multiple  Hypothesis 
Testing  in  Space  Domain  Awareness  with  the  Space  Surveillance  Telescoped  Applied 
Optics.  Accepted  April  2016 

Conference  Papers  and  Presentations: 

•  Based  on  Chapter  5: 

T.  Hardy  and  S.  Cain.  “ Characterizing  point  spread  function  (PSF)  fluctuations 
to  improve  resident  space  object  detection  (RSO)d  in  SPIE  Defense-i-  Security, 
International  Society  for  Optics  and  Photonics,  2015. 

•  Based  on  Chapter  4: 

T.  Hardy  and  S.  Cain.  “ Investigating  prior  probabilities  in  a  multiple  hypothesis 
test  for  use  in  space  domain  awarenessd  in  SPIE  Defense-r  Security,  International 
Society  for  Optics  and  Photonics,  2016. 
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